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ABSTRACT 

We investigate a model in which galactic nuclei form via the coalescence of pre-existing stellar systems 
containing supermassive black holes. Merger simulations are carried out using A-body algorithms that 
can follow the formation and decay of a black-hole binary and its effect on the surrounding stars down 
to sub-parsec scales. Our initial stellar systems have steep central density cusps similar to those in 
low-luminosity elliptical galaxies. Immediately following the merger, the density profile of the remnant 
is homologous with the initial density profile and the steep nuclear cusp is preserved. However the 
formation of a black-hole binary transfers energy to the stars and lowers the central density; continued 
decay of the binary creates ap~r"' density cusp similar to those observed in bright elliptical galaxies, 
with a break radius that extends well beyond the sphere of gravitational influence of the black holes. 
Our simulations are the first to successfully produce shallow power-law cusps from mergers of galaxies 
with steep cusps, and our results support a picture in which the observed dependence of nuclear cusp 
slope on galaxy luminosity is a consequence of galaxy interactions. We discuss the implications of our 
results for the survivability of dark-matter cusps. 

We follow the decay of the black hole binary over a factor of ~ 20 in separation after formation of 
a hard binary, considerably farther than in previous simulations. We see almost no dependence of the 
binary's decay rate on number of particles in the simulation, contrary to earlier studies in which a lower 
initial density of stars led to a more rapid depletion of the binary's loss cone. We nevertheless argue that 
the decay of a black hole binary in a real galaxy would be expected to stall at separations of 0.01 — I pc 
unless some additional mechanism is able to extract energy from the binary. 



1. INTRODUCTION 



Galactic nucleQ are regions of high stellar density at 
the centers of g alaxies. Early studies of the evolution of 
galac t ic nuclei ([Spitzcr fc Saslaw. 1966 ; Spitzcr & Stone 



1967; Colgate, 1967; Sanders, 1970) emphasized stellar en- 
counters and collisions as the dominant physical processes. 
In these models, the density of a compact stellar system 



about 1990 lacked the resolution to determine whether t r 
was shorter than a Hubble time on scales smaller than ~ 1 
pc in galactic nuclei. We now know that stellar densities 
increase approximately as power laws in galactic nuclei, 
p ^ r~ 7 , down to the smallest radii that can be resolved, 
or roughly 10 _1 pc in the nearest galaxies. For instance, 
the Local Group galaxies M31, M32 and the Milky Way all 



gradu ally increases as energetic stars arc scattered into 



have nuclear density cusps with 7 = 1.5±0.5 (Lauer et al 



elong ated orbito via two body oncountoro. The incroaoo in 
density leads to a higher rate of physical collisions between 
stars, liberating gas that falls to the center of the system 
and c ondenses into new stars w hich undergo further colli- 
sions. Begelman fc Rees (1978 ) argued that the evolution 
of a dense nucleus would lead inevitably to the formation 
of a massive black hole (BH) at the center, either by run- 
away stellar mergers or by creation of a massive gas cloud 
which collapses. Subsequent studies ( Duncan fc Shapi ro 



1998). Furthermore the evidence for supermassive BHs is 



compelling in these galaxies. Within the BH's sphere of 



influence r gr , where 
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the stellar velocity dispersion rises as 
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(0) then implies t r oc r 7 3 / 2 w r°, nearly independent 
of radius. The relaxation time based on observations at 



0.1" exceeds 10 11 yr in almost all galaxies (Faber et al 



eluded "seed" riris which grow via accretion ot stars or gas 
liberated by stellar collisions. 

The fundamental time scale in these models is the re- 
laxation time determined by the stars, or 



1997); this angular size corresponds roughly to a radius 
for nearby galaxies ( Merritt fc Fcrrarcsc, 2001b| ) , from 



0.34 



flSpiti 



G 2 m*phiA 



(1) 



tzer & Mart, lyYl), where cr* is the stellar 1L» velocity 
dispersion, to* and p are the stellar mass and mass density, 
and In A is the Coulomb logarithm. Observations before 

1 We use the term "nucleus" to refer generically to the central 
parts of galaxies. The term is sometimes used more restri r.tively to 



' gr 

which it follows that t r is likely to exceed a Hubble time 
at smaller radii as well. ( The pointlikc nucleus of M33 is 
probably an exception; see Kormcndy & McClure (1993).) 

Physical collisions between stars occur on a timescale 
that is longer than t r by roughly a factor (In A )8 2 /(l 
8) 10 where O, the "Safronov number" (safronov, 
1960), is of order a few for stars in a galactic nucleus 



refer to pointlike nuclei, e.g. Kormendy & McClure (1993) 



Thus neither elastic nor inelastic gravitational encounters 
are likely to be of dominant importance in determining the 
structure of nuclei containing supermassive BHs. 

Nevertheless the properties of galactic nuclei do vary in 
systematic ways with the properties of their host galaxies 
( Kormendy, 1985 ; Lauer, 1985 ; Faber et al., 1997 ) and one 



2 



would like to understand this. Recent discussions of the served in the brightest elliptical g alaxies flMerritt fc Frid- 



forma tion and evolution of galactic nuclei have begun from 
the assumption that supermassive BHs were created dur- 
ing the quasar epoch and have been present ever since with 
roughly their current masses. Another element missing 
from the earlier studies was galactic mergers. Mergers are 
complex phenomena, but an almost certain consequence of 
a merger is the infall of the prog enitor galaxy's BHs into 
the n ucleus of the merged system ( Begelman, Blandford & 
Rees, 1980| ). An infalling BH would be expected to carry 
with it a mass in stars of order its own mass, and decay 
of the BHs' orbits would inject a substantial amount of 
energy into the stars, enough to determine the structure 
of t he remnant nucleus out to a radius of several times 
(Ebisuzaki, Makino & Okumura, 1991). In this pic- 



' gr 

ture, the structure and kinematics of galactic nuclei are 
fossil relics of the merger histories of galaxies and of the 
interaction between stars and supermassive binary BHs. 

The present paper is a numerical study of this forma- 
tion model. We simulate the merger of two galaxies, each 
of which contains a central point mass representing a su- 
permassive BH. Our study is unique in that it follows the 
details of the merger from its earliest stages, when the two 
galaxies are distinct, to its late stages, when the BHs have 
formed a hard binary and the binary has decayed via en- 
ergy exchange with stars to a separation much less than 
one parsec. No existing iV-body code can efficiently fol- 
low the evolution over such a wide range of scales; hence 
we break the calculation into two parts, before and after 
formation of the BH binary, and use different algorithms 
for each (§2). 

We also include for the first time initial models which 
are self-consistent realizations of galaxies with steep cen- 



tral density cusps, p > 



This choice is motivated by a 



number of lines of argument which suggest that the density 
of stars around a supermassive BH should be a steep power 
law. Random gravitational encounters between stars lead, 
over two-body relaxation times, to density profiles of the 
2 23 in the absence of a black hole flCohn 



for 



198(f and p ~ r in the presence of a BH ( pahcal 
& Wc lf, 197g ). It was argued above that relaxation pro- 



cesses are probably of secondary importance in most nu- 
clei, but even the slow growth of a BH in a pre-existing, 
collisionless nucleus produces a density pro file of the form 
p ~ r~ 7 , r < r or with 1.5 <; 7 < 2.5 (|Peebles, 1972|; 
Young, 1980|; Quinlan, Hcrnquist fc Sigurdsson, 1995 ; |vai]| 



der IV arcl, 19991 ). Low-luminosity elliptical galaxies and 
the bulges of spiral gal axies are observed to have steep 



Lauer et 



cusps, 1.5 <; 7 <; 2.5 ( Fcrrarese et al., 1994 
al., 1£ 95); these galaxies are the least likely to have been 
strongly affected by mergers and hence their density pro- 
files may reflect the structure of all nuclei at early times. 
Finally, hierarchical growth of structure in the universe 
generically produc es systems with steep centr a l density 



cusps. 



-1.5 



( Dubinski fc Carlberg. 1991 ; Navarro, 



Frenk |fc White, 199lpVloore et a i, ; 1998|) , although sim- 
ulations do not yet have sufficient resol ution to make pre- 
dictions on parsec or sub-parsec scales ( Moore, 2001 ). 

A major success of our study is the demonstration (§|]]4|) 
that the merger of two galaxies with steep, power-law den- 
sity cusps can produce a galaxy with a shallow power-law 
cusp. Shallow cusps (also called "cuspy cores") are ob- 



man, 199E; Gebhardt et al., 1996), and, while their origin 



has tentatively been associated with the binary BH mode l 
( |Ebisuzaki, Makino fc Okumura, 199lfc [Faber et al., 1997| ), 
no previous simulation had the resolution necessary to fol- 
low the coalescence of initially steep cusps. Our results 
support a model in which the observed dependence of nu- 
clear cusp slope on galaxy luminosity is a consequence of 
galaxy interactions (§||). 

We also discuss in detail (§Q) the decay of the BH binary; 
we follow that decay over a factor ~ 20 in semimajor axis 
after formation of a bound pair, considerably farther than 
in earlier simulations. An important question is the de- 
pendence of the binary hardening rate on TV, the number 
of particles in the simulation. Earlier studies had noted 
a decreasing decay rate with increasing N, implying that 
the decay in real galaxies, where N is very large, might 
be slow. We do not observe an appreciable N dependence 
in our simulations; we discuss the likely reasons for this in 
§||, but argue nevertheless that the decay would probably 
stall in real galaxies, at separations of 0.01 — 1 pc, unless 
some additional mechanism is effective at extracting en- 
ergy from the binary. We are therefore led to predict (§^) 
that some galaxies contain uncoalesced BH binaries at the 
current epoch. 

We also present the morphological and kinematical struc- 
ture of the merged galaxy on sub-parsec scales (§||) and 
discuss some observational signatures associated with our 
formation picture. 

2. METHOD 

In this section we describe the initial conditions and 
algorithms used in our simulations and compare them with 
those of other authors. 

Initial conditions consisted of twin, spherical stellar sys- 
tems following the density law 



p(r) = 



M 



r 

1 + — 

7-0 



(3) 



( |Jaffc, 1983 ; Dehncn, 1992), where tq is the half-mass ra- 
dius and M the total mass. To each of the models was 
added a central point of mass M, = 0.01 M represent- 
ing the supermassive BH. The ratio M,/M in our models 
is somewhat greater than the mean ratio of B H mass to 
galaxy mass in ob served galaxies, ~ 1.2x 10 -3 (Merritt & 



Ferrarese, 20011:); however the radius of influence of our 
BHs in our merged galaxies, r gr — GM,/a1 ~ 0.02, is 
much smaller than ro which allows us to ignore the large- 
scale stellar distribution when scaling our models to the 
nuclei of real galaxies; see §||. Velocities of the stars were 
generated from the uniq ue, isotropic distribution function 



( Tremaine et al., 1994 ) that reproduces the density law 
(|) in the combined potential of the stars and the cen- 
tral point mass. Thus our models are initially in a state 
of detailed equilibrium. The values chosen for the model 
parameters are listed in Table |l|; Newton's constant is set 
to unity.Q The galaxies were set at time t = in an el- 
liptic mutual orbit of semimajor axis ac = 2 and velocity 
vq — 0.1425 equal to half of the circular orbit velocity. 



1 These parameters agree with the Heggie & Mathieu (1986) "stan- 
dard" units. 
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Table 1 
Model Parameters 



Parameter 


Symbol 


Value 


Mass of Galaxy 


M 


1 


Mass of Black Hole 


Mm 


0.01 


Half-Mass Radius 


ro 


1 


Total Energy 


E 


-0.25 


Number of Stars 


N 


131,072 


Stellar Mass 


m* 


7.63 x 10- 6 



(The index "G" labels the binary composed of two galax- 
ies as distinct from a binary composed of two BHs.) Since 
the relative orbit rapidly circularizes, we do not expect 
our results to be strongly affected by the choice of orbital 
initial conditions. 

Our goal was to follow the evolution of this binary sys- 
tem from its earliest stages, when the galaxies were dis- 
tinct, through the formation of a bound BH pair, until 
the gradual exchange of energy between BHs and stars 
had caused the BH binary to shrink to sub-parsec separa- 
tions. No single TV-body code currently in existence can 
deal efficiently and accurately with evolution over such a 
broad range of scales; particularly demanding is the treat- 
ment of BH-BH and BH-star interactions, which require 
an algorithm that can accurately integrate the equations 
of motion of point masses without softening. The clos- 
est approximation to such a code is SCFBDY developed 
by G. Quinlan at Rutgers University and u s ed in two 
publi shed stud ies (Quinlan & Hcrnquist, 1997; Merritt & 
Quinlm, 1998). SCFBDY combined elements of Aarseth's 
NBODY series of codes, including regularization of the BH- 
BH in teraction, with a low-resolution, mean-field potential 



solver for computing the force field due to the stars. How- 
ever SCFBDY is only suited to systems with a single den- 
sity center and a high degree of symmetry, ruling out its 
application to mergers. SCFBDY also uses softened gravity 
for co mputing the BH-star interactions, an approximation 



that affects the accuracy of the critical interactions leading 
to the decay of the BH binary. This code is also not avail- 
able in a form that runs on parallel architectures, limiting 
the number of particles that can be used. 

We chose to break the problem into two parts, using a 
tree code for the early stages of the merger (roughly un- 
til the formation of a BH binary), and a high-precision, 
direct-summation code for the later stages. Ideally, one 
would use a high-precision code right from the start, to 
handle the steep force gradients produced by the BHs and 
the other stars in the cusps. However we were willing 
to accept some inaccuracy in the integrations during the 
early stages of the merger if in so doing we could treat a 
larger N; our only prerequisite was that motions of stars 
on scales larger than the separation of the hard BH bi- 
nary should be faithfully tracked during the early stages 
of the merger. This reasoning motiva ted our choice of the 
recen tly- relea sed tree code GADGET (3pringcl, Yoshida & 
White, 2001) as the integrator for the early stages of the 
merger. Features of GADGET that are relevant here in- 
clude domain decomposition of the particle data set, map- 
ping of particles onto the classic octal tree structure that 



respects hierarchical clustering of particles, quadrupolar 
expansion of force moments for spatially separated nodes, 
and individual and adaptive time steps for all particles. 
Force integration is controlled through the parameters 
(h,rj) that denote, respectively, the gravitational soften- 
ing length and the time step accuracy factor, At — rj/\a\. 
Special care was taken to identify parameter values lead- 
ing to optimum accuracy and efficiency on the Rutgers Sun 
HPC-10000 and the SDSC Cray T3E systems where the 
GADGET runs were produced. The softening length was 
chosen to be h — 0.001, smaller than both the BH gravi- 
tational radius r gr = GMnjo^ as 0.01 and the separation 
corresponding to a hard binary, = GMi 2 /8cr 2 fts 0.0025. 
We monitored the density profiles and Lagrangian radii 
of the stellar cusps as diagnostics sensitive to corrup- 
tion of the bulk stellar distribution due to unacceptable 
levels of softening. Runs entering final selection exhib- 
ited no cumulative distortions on scales larger than ~ h. 
The total particle number used in the production run was 
N = 2 18 = 262, 144. 

The late stages of the evolution were int egrated us- 
ing the direct-sum mation code NBODY6++ (3purzem & 



Baumgardt, 1999). Conceived for the study of relaxation 
phenomena in globular clusters, NBODY6++ and its se- 
rial progenitor NBODY6 are th e last and mos t complex 
codes in the NBODYx series of Aarseth (1999|), employ- 
ing th e fourth-order Hermite scheme (Makino & Aarseth 



1992) as their primary integrator. The codes were writ- 



ten to facilitate the exact integration (no softening) of a 
large number of bodies with approximately equal masses. 
NBODY6++ gives particles adaptive block- individual time 
steps At n tx 2 _n that are short for particles in dense re- 
gions and as much as 10 2 time s longer for particles at 
the outskirts of the system. The Ahmad & Cohen (1973) 
scheme is used to select a subset of neighbors whose forces 
on the test particle are extrapolated at higher time reso- 
lution than those of the non-neighbors. Near encounters 
are treated using the two-body KS regularization scheme 
(Kustaanheimo & Sticfel, 1965) and its generalizations to 
systems with a few bodies, including the triple, quad and 
chain regularizations. Chain regularization is a system- 
atic procedure for serializing the pairwise KS regulariza- 
tion in a group of not more than six bodies in close ap- 
proach. NBODY6++ was developed to run on low-latency 
distributed memory architectures such as the Cray T3E. 
There is no domain decomposition: every processing node 
contains an identical copy of the whole dataset. Only 
the do-loops are broken into parallel sections; after every 
force calculation, an all-to-all broadcast scheme updates 
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Table 2 
NB0DY6++ Initial Conditions 



Label 


Reduction in Particle Number 


Truncation in Energy 


N 


M./m* 


GADGET 


lx 


lx 


262, 144 


1,311 


A2 


2x 


4x 


32, 768 


655 


A4 


4x 


4x 


16, 384 


328 


A8 


8X 


4x 


8, 192 


164 


B2 


2X 


4x 


32, 768 


1,311 


B4 


4x 


4x 


16, 384 


655 


B8 


8x 


4x 


8,192 


328 



the particle sets. The treatment of binaries (regulariza- 
tion etc.) is not parallel. In the simulations conducted on 
the 272-node Cray T3E at the SDSC and the 64-node Sun 
HPC-10000 at Rutgers, NB0DY6++ scaled well with the 
number of processors when the spread in time steps was 
moderate. When a few particles had much shorter time 
steps than others, the scaling was poor and serial integra- 
tion on 666 MHz Alpha chips was found to be preferable. 

One goal of our study was to identify any iV-dependent 
features of the evolution, since real galaxies have much 
larger numbers of stars than accessible to direct simula- 
tion. To isolate the effects of varying N from other depen- 
dence (and also because NB0DY6++ can not deal with 
particle numbers as large as 10 5 on currently available 
machines), we drew various random particle sets from the 
N = 2 18 GADGET integration and used these as initial 
conditions for the NB0DY6++ runs. Half of the particles 
were iteratively removed and the masses of the remaining 
particles (except for the BHs) were doubled. We chose the 
time to = 10.6 to select our reduced data sets; the separa- 
tion of the BHs at this time is 0.072, substantially greater 
than their separation (~ 0.0025) at the time th when a 
hard binary forms. 

In addition, to increase the effective resolution near the 
center, we sorted the data set by energy Ei — m*[uf/2 + 
U(Yi)\ and removed the upper 3/4 of all stars. While the 
new data sets were statistically distinct from the original, 
we expect to hnd unchanged dynamics in the cusps where 
the fractional perturbation from the removed stars is neg- 
ligibl e. Combination of these two techniques led to a set of 



ular cluster simulation were found to be weaknesses when 
the code is applied to systems with massive BHs. Chain 
regularization becomes impractical when some of the bod- 
ies are much more massive than others. In a model galaxy 
with N = 10 6 stars and a dense p ~ r~ 2 stellar cusp har- 
boring a M./m* = 10 3 central BH, there are ~ 10 3 bodies 
whose orbits are largely determined by the BH's potential. 
Although these stars satisfy the requirements for KS regu- 
larization (tiny mutual separations and short time steps), 
only the forces between the BH and the stars are signif- 
icant; star-star encounters are energetically unimportant. 
However NBODY6++ is incapable of making this distinc- 
tion. It will either try to generate a KS chain containing 
10 3 bodies, which would defeat the purpose of regulariza- 
tion altogether and is beyond the present design; or, if one 
"turns off" the chain regularization, NBODY6++ will re- 
sort to a pairwise KS regularization whenever two stars 
come close, even if the motion of both stars is mostly in 
response to the force from a BH. In order to make the 
integrations go efficiently, it was sometimes necessary to 
remove "by hand" a few particles that were in tight orbits 
around one of the BHs. 

Our simulations uniquely incorporate three features. 1. 
The galactic merger leading to the formation of the bi- 
nary BH is carried out from a state in which both galaxies 
are reasonably isolated. 2. Both galaxies initially contain 
faithful realizations of steep stellar density cusps. 3. The 
BHs are present in the cusps at the outset in a kinemati- 
cally consistent manner. Earlier studies have in corporated 
some of these features but never all of them. Makino & 



formula "An" where n is the fractional reduction (Table I 

In order to verify the accuracy of the tree code inte- 
grations, we continued the GADGET run for several dy- 
namical times after to and compared the results with the 
NBODY6++ integrations. Coincidence was found to per- 
sist until the separation 7*12 between the BHs was ~ h, at 
which point the BHs start to "see" each other's spurious 
finite extent in the softened GADGET integrations (Figure 
[l]). From that point on, the binary separation saturates 
with GADGET but continues to decrease with NBODY6++ 
. We conclude that GADGET faithfully reproduces the dy- 
namics of merging stellar cusps and BH binaries in regimes 
where the binary is soft. 

NBODY6++ would seem to contain all of the machinery 
necessary for efficiently and accurately handling star-star, 
BH-BH and BH-star interactions. In fact, some of the 
greatest strengths of NBODY6++ in the context of glob- 



initial conditionsforNBO^ Ebisuzaki (1996) reported mergers of King models with 



N = 16,384 particles and central BHs with M. > 1/64. 
They conducted repeated mergers by recycling the merger 
products into initial conditions for successive mergers, af- 
ter replacing the pair of BHs by a center-of-mass particle 
and reducing the number of particles by half. Makino 
& Ebisuzaki's choice of King models as initial conditions 
made their galaxies poor representations of real stellar 
spheroids which always contain power-law density cusps; 
nor could they test the hypothesis that weak cusps are 
generated by t he interaction of BHs and surrounding stars. 



Makino (1997) studied the evolution of a BH binary pro- 
duced by a similar sequence of mergers, also using King- 
model initial conditions, but with a much larger maximum 
number of particles, N = 262, 144. While the number 
of particles in Makino's and our simulations is effectively 
the same, judging from the density profiles in Figure 4 of 
his paper, Makino's initial conditions appear to be ~ 10 2 
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Fig. 1. — (a) Evolution of the separation between the BHs. Time is measured from the start of the simulation in units such that the crossing 
time in a single galaxy is ~ 2.2. The BHs form a hard binary at t = 11.0. Squares are from the 2.62 X 10 5 -particle integration with the 

tree code GADGET ; the binary separation saturates at roughly the softening length h, marked by the arrow, in this simulation. Solid line is 
from the NBODY6H — h integration A2 with M./m* = 655. NBODY6H — h is able to follow the decay of the binary to arbitrarily small scales, 
(b) Stellar density as a function of time of stars separated by distance r < 0.04 from either of the BHs. (c) Stellar velocity dispersion within 
0.01 < r < 0.04 around each black hole. 



times less dense than ours inside of the binary's sphere of 
influence. Quintan & Hcrnquist (1997) studied the evo- 
lution of a BH binary inside cuspy models with p ~ 



and p 



and a wide range of BH masses and particle 



numbers, N < 2 x 10 5 . As discussed above, Quinlan's code 
was unable to simulate an actual merger due to the limita- 
tions of its treatment of the mean field. All of the detailed 
results in their paper were derived from initial conditions 
consisting of a single galaxy into which two "naked" BHs 
were dropped from starting points located diametrically 
apart at the half-mass radius. This configuration is likely 
to produce substantial evolution of the cusp before the for- 
mation of the binary as the infalling BHs heat the stars; 
this can in fact be seen in their Figure 1, a plo t of La- 
grang ian radii over time . As in the sim ulations of Makinc 
& Ebi|suzaki (1996|) and[Makino (1997[), the initial density 



of sta rs around the BH binary in the simulation of Quin 
lan & Hcrnquist (1997 ) was much lower than in ours. This 



difference will t urn out to be consequential (§4). 

Barnes (lyyy) presented simulations of mergers of iden- 
tical galaxies with power-law cusps and no BHs. He 
used N — 65, 536 bodies and a fixed time step, leap- 
frog scheme. Barnes showed that, outside of the softening 
length, power-law cusps as steep as 7 = 2 are preserved 
by the merger. We find an analogous result (Q. 

Our simulations, being purely dynamical, preclude any 
non-dynamical processes such as gas-driven dissipation 
that might act to accelerate the binary's decay. Decay 
might also be enhanced by dynamical processes that we do 
not include, e.g. the passage of a star cluster, gas cloud or 
third supermassive BH through the nucleus. How different 
is the signature on the stellar distribution of a binary that 
coalesces immediately after its formation? We address this 
question by shadowing each simulation introduced above 
with another where the BH binary is replaced by a single 
BH of mass 2M. at time t = t + 0.4 = i/,. This yielded 
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Fig. 2. — Projected density contours for the run A2 (NB0DY6H — h ) with Af./m* = 655. The orbital motion of the BHs (positions indicated 
by filled circles) is clockwise in the plane of the figure. First row: t = 10.67, 10.8, 10.98, 10.91. The two cusps spiral-in under the influence of 
dynamical friction. Second row: t = 10.94, 10.95, 10.96, 10.97. The two cusps merge into one. The final density profile is similar to that of the 
initial stellar systems. Third row: t = 10.98, 11.0, 11.1, 11.6. The density of the newly-formed cusp drops rapidly as the BH binary transfers 
energy to the stars. Fourth row: t = 12.6, 13.6, 16.6, 18.6. Density continues to drop as the BH binary ejects stars. 



a second set of initial conditions labelled as "Bn" in Table 



3. CUSP COALESCENCE 

We divide the evolution into two regimes, before and af- 
ter the formation of a hard BH binary, and discuss the first 
regime in this section. The two regimes correspond ap- 
proximately, but not exactly, to the intervals before and af- 
ter the start of the NBODY6++ integrations; as discussed 
above, these integrations were begun when the BHs were 
still a few softening lengths apart and had not yet formed 
a tig htly bound pair . 

As Quinlan (1996) notes, there are many ways to define 
a "hard" binary. The standard definition is energy based: 
a binary is hard if its binding energy exceeds the typical 
particle kinetic energy, \Eb\ 3m»er 2 /2. This definition 
is inapplicable to the case of massive BH binaries since the 
binary-to-stellar mass ratio scales with N while a is inde- 
pendent of N; according to this definition, a very massive 
BH binary would always be hard if bound and soft oth- 



erwise. The famous law of Hcggic (1975|) , asserting that 
hard binaries evolve toward even harder states, suggests a 
second definition of hardness. While a viable distinguish- 
ing criterion for hard binaries in star clusters, Heggie's law 
fails to capture the transition between two different pro- 
cesses — dynamical friction and mass ejection — that both 
tend to drive a massive binary to an ever-harder state in 
our simulations. 



We therefore followed the suggestion of Hills (1983) and 
Quinlan (1996) and defined a "hard" binary in terms of 
its orbital velocity. The orbital velocity of each BH in 
a circular-orbit binary is v 2 c = GM^/Aa with Myz the 
combined mass of the two BHs and a their separation. We 
defined the critical separation at which a binary becomes 
"hard" as 

GM 12 

ah = ^f (4) 

corresponding to v c = y/2o~*. In model units, ah ~ 
2.5 x 10~ 3 and the binary separation first falls below ah 
at t = t h w 11.0 (Figure p. "Subsonic" massive binaries, 
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a> ah, harden by dynamical friction acting on each BH 
(and its associated cluster) individually; a "supersonic," 
or hard, binary behaves like a structureless point mass un- 
der the action of dynamical friction but can capture stars 
and eject them at much higher velocity, thereby increas- 
ing its hardness. Quinlan (1996| ) notes that this definition 
of hardness is roughly equivalent to the statement that a 
hard binary hardens at a constant rate. We found this to 
be true (§4); the definition (|) is also a natural one in the 
sense that the character of the binary's evolution, and the 
evolution of its surrounding stellar cluster, were found to 
undergo qualitative changes when a dropped below ~ a/,, 
corresponding to the onset of mass ejection by the binary. 

The merging of the stellar cusps for run A2 is illustrated 
in Figure g. This figure makes manifest that the BHs re- 
main closely associated with their initial stellar cusps dur- 
ing every stage of the merger, up to and including the point 
when the two cusps merge into one at t « t/,. A conse- 
quence is that dynamical friction brings the BHs together 
much more rapidly than if they were "naked, " since their 
effective mass is the mass of the cluster of stars bound to 
them. 

We can check this assertion by comparing the orbital 
decay rates for an isolated BH and that embedded in a 
stellar cusp. The decay rate for an is olated BH on an ap 



proximately circular orbit is given by ( Binncy & Trcmainc 
19870 



da 
~dt 



erf(l)-erf'(l)GM. 



-0.302 



V2 
GM. 



In A 



o~*a 



In A 



(5) 



where a denotes separation between two BHs. With 
M. = 0.01, a = 0.1 and er* »s 2 -1 / 2 , the decay rate is 
estimated at da/dt ss — 0.043 In A. Note that for a > 0.1 
the formula yields even lower rate of decay. As for the 
Coulomb logarithm, it can be written in terms of the ra- 
tio of the maximum and the minimum impact parameter 
In A = \n(jj max i /pmin) where it is a standard choice to se- 
lect the gravitational radius of the isolated "test particle" 
for the latter, p m i n = GM,/cr%. In lack of a canonical 
choice for p ma x, we equate it to the orbital radius, which 
implies InA = ln5.0 « 1.6 and thus da/dt ~ —0.07 (but 
see also Appe/ndix A). The predicted decay rate is a fac- 
tor of ~ 6 smaller than the rate of in-spiral da/dt —0.43 
we measured in the GADGET run in the interval 9.0 < 
t < 11.0 preceding the formation of hard binary, in this 
interval a(t) = 0.78 - 0.43 X (t - 9.0) is a good fit. 

We compare the isolated particle estimate with an esti- 
mate of how rapidly dynamical friction would act to bring 



together two overlapping spheres with p 



density 



profiles. Let the spheres each have mass A4 and density 
p = cr 2 /2irGr 2 inside a radius a. When the separation a 
between their centers is much larger than the BH radius 
of influence, or equivalcntly M.(a) S> M., we can ignore 
the BHs. Then At (a) — a 2 a/G and the circular velocity 
is v 2 (a) = a 2 /2. The dynamical friction force acting on 
one of the spheres is given by 



(A«n) = 
F(v) = 



±TrG 2 Mp\nAF{v) 



erf(x) 



ccerf '(x) 



( Chandrasokhar. 1943| ) , where InA is the Coulomb loga- 
rithm and x = v / \/2cr* . Setting v — v c — er*/\/2 gives 
F = 0.0811. Taking for p the density of either sphere at a 
distance a from its center, we find 



(Awn 



0.324cr? InA 



(7) 



Equating the torque produced by this acceleration with 
the change of orbital angular momentum J, and writing 



dJ dJ da 
dt da dt 



<7 J a da 
V^Gdi' 



gives 



da 



-0.23ct» InA. 



(8) 



(9) 



Just outside of the BHs' sphere of influence, cr* « 1 (Figure 
0b). For InA we again take ~ 1.0, giving da/dt « 0.24. 
This result still falls short of the observed value (da/dt « 
-0.43) by a factor of - 2. 

A detailed integration (Appendix A) taking into account 
the shape and the finite extent of both spheres yields the 
rate (A«||) w — 1.50a 2 /a (note the absence of InA), there- 
fore da/dt sa —1.06. This result, however, is sensitive to 
the choice of tidal radius outside of which the spheres are 
indistinguishable; estimation of this radius is difficult in 
part due the large orbital eccentricity of the galaxies in 
the simulation. 

Remarkably, the density structure of the merged galaxy 
just after formation of a hard binary is essentially identical 
to that of the initial stellar systems at radii r p> a. This is 
illustrated in Figure |[ Homology following a merger was 
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Fig. 3. — Radial density profiles of the pre-merger galaxies (thin 
curve) and at time t ss t h = 10.96 (thick curve) when the binary 
separation equals ri2 = 6.5 X 10 — 3 . The pre-merger density was 
multiplied by the factor (Mr^ ^)new /(Mr^ 3 ) old as 0.53 to bring 
it to the scale of the post-merger galaxy. The two galaxies have 
merged into a single galaxy that is nearly homologous with the initial 
galaxies on scales r > r\%/1. Shortly after this time, the central 
density drops as the binary heats the core. 
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Fig. 4. — Lagrangian radii around each BH in the first time unit 
of the NB0DY6H — h run A2. From bottom to top, the radii enclose 
1CT 4 , 10" 3 ' 5 , lO- 3 , 1CT 2 ' 5 , 10" 2 , lO" 1 ' 5 and 1CT 1 in units of the 
mass of one galaxy before the merger. 
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found also by Barnes (1999) in spherical galaxies without 
BHs. In our simulations, however, the homology is short- 
lived. The formation of a hard BH binary at t « th is 
followed by a sudden drop in the stellar density within the 
binary's gravitational sphere of influence, r < 0.01. This 
is clearly seen in Figure |], a plot of Lagrangian radii, and 
also in the density contour plot of Figure ^. (The drop in 
density is not so apparent in Figure |]b because the density 
plotted there is an average over a radius of 0.04, and little 
net change in density occurs within this radius - see Figure 
) In effect, the steep cusp that was present immediately 
after formation of the BH binary is destroyed in little more 
than the local crossing time. 

What is responsible for the rapid destruction of the 
cusp? Two possible, and closely related, mechanisms are 
deposition of energy into the stars by dynamical friction 
acting on the BHs individually; and ejection of stars that 
exchange energy with the BH binary. Neither process is 
well defined in this regime where the BH binary is nei- 
ther very hard nor very soft. Nevertheless we can write 
approximate expressions for the rate at which energy is 
transferred to the stars by the two mechanisms, by assum- 
ing either that the BHs are moving independently of each 
other, or as members of a tight binary; of course neither 
assumption is strictly satisfied. 

In the first case, dynamical friction would extract energy 
from the two BHs at a rate 



(AE) = 2Mv(Av\\ 



SnG 2 M 2 p\nkF(v) 



ejection, produces energy at the rate 
dE G 2 M 2 pH 



dt 



2<t* 



(12) 



( [Hills, 19 831 
where H 



Quinlan, 1996), 



Mikkola fc Valtoncn 1992] ^ 
is the dimensionless hardening rate; H ~ 15 m 
the limit of a very hard, equal-mass binary and drops to 
~ 10 for a binary with a = ah- Thus, both mechanisms 
predict an energy deposition rate that can be written as 



dE 
~dt 



CG 2 M 2 pa, 



(13) 



with C w 5; we have taken In A sa 0.5 (Appendix A). 

The energy of the binary when a = ah is E = — 2Ma 2 , so 
the characteristic time over which either process extracts 
energy is ~ 2C~ 1 a^/G 2 Mp. Computing p using the mass 
within r = 0.01 at t = th gives an energy extraction time 
scale of ~ 0.2 in model units. This is quite comparable to 
the time associated with the jump in Lagrangian radii of 
Figure ||. We note also that the energy extracted from the 
binary in this time, E « 2Ma 2 0.02, is comparable to 
the energy in stars within a radius of ~ 0.01. This too is 
consistent with the changes in Lagrangian radii shown in 
Figure ^. 

We conclude that the sudden disruption of the steep 
cusp is attributable to transfer of energy from the BHs into 
the surrounding stars as the BHs form a hard binary. We 
emphasize again the rapidity of this process, which earlier 
analyses have overlooked. If a galaxy's cusp is to avoid 
this fate, some mechanism must extract energy from the 



Y 



(10) 



-0.1 




(cf. equation H) . Setting v — our definition for the 

onset of a hard binary, gives F = 0.43~l/2 and 



(AE) w -2V2irG 2 M 2 po-- 1 lrLA. 



(11) 



The alternate mechanism, hardening of the binary by mass 



Fig. 5. — Isophotes in the run A2 in three projections. Black dots 
show the location of the BH binary; separation between the BHs at 
this time is a ~ 1.5 X 10~ 4 . 100 snapshots of the nucleus were 
superposed in the interval 18.1 < t < 19.1. The merger remnant 
is approximately axially symmetric with an edge-on ellipticity of 
6 w 0.25. 
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-0.2 -0.1 0.1 0.2 

Fig. 6. — 2D kinematics of the merging cusps for the run A2. View is in the plane of the merger from a direction perpendicular to the 
line connecting the two BHs. Left panels show the mean line-of-sight velocity; blue contours indicate approaching stars. Right panels are 
the line-of-sight velocity dispersion. In all panels the contours are separated by 0.038. (a) t = 10.66; (b) t = 10.82; (c) t = 10.96. The BHs 
remain centered on the velocity dispersion peaks but move inward with respect to the peak of the rotational velocity. 
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Fig. 7. — Spatial density profiles (left column) and projected density profiles (right column). In each panel, starting at time t = 11.0, the 
profiles are recorded from top to bottom at intervals At = 1.0, each of them an average obtained by superposing 100 snapshots. Top row, 
run A8; middle row, run A4; bottom row, run A2. Thick lines represent the run B2 with one BH; the merger remnant would have this profile 
if the BHs coalesced at t = t h . Dashed lines are profiles of the original galaxies multiplied by an arbitrary factor. The merger remnant has a 
p ~ r _1 density cusp which projects to a core profile with continuous curvature. 
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binary BH on a time comparable to the local dynamical 
time, before it is able to exchange energy with the stars. 

The density near the BHs continues to drop at later 
times although more slowly, as the BH binary gradually 
decays. We discuss this process in more detail below. 

The large-scale kincmatical evolution of the merger is il- 
lustrated in Figure |[ The general character of these plots 
is similar to what is seen in si mulated mergers of cqua l- 
mass galaxies without BHs (e.g. Bendo fc Barnes (2000D ), 
with the highest mean rotational velocities in the plane 
of the merger and a roughly cylindrical rotation pattern 
elsewhere. The peaks in the rotational velocity initially 
correspond with the locations of the BHs, but dynamical 
friction causes the BHs to move inward from these peaks 
at a time t w 10.75. Velocity dispersions, on the other 
hand, remain peaked on the BHs at all times consistent 
with the fact that the BHs remain centered on their cusps 
(Figure ||) . At the time of formation of the hard BH bi- 
nary, t — th « 11.0, the merger remnant is mildly rotating 
with a peak line-of-sight rotational velocity of v/a* « 0.58 
at a distance ~ 0.04 from the center (using w 0.6 at 
distance 0.2), although bulk fluctuations at 20% level in 
the rotation field persist on scales r <; 0.1. The ellipticity 
of isophotes shown in Figure | is e « 0.25 w hich falls on 
the isotropic oblate rotator relation v/er* = ye/Jl — e). 

Evolution of the stellar density profiles is shown in Fig- 
ure 0. Profiles were computed from the A-body positions 
using a n onparametric kernel routine based on the algo- 
rithms in Merritt fc Tremblay (1994 ); details are given 
in Appendix B. Each profile is an average over several 
snapshots; they are separated by At = 1.0 starting at 
t = th = 11.0. The evolution of p(r) becomes more regu- 
lar as N is increased, due probably to the smaller random 
motion of the binary for larger N. Considerable evolu- 
tion occurs in the interval th — 1 <J t <^ th + 1 when the 
BHs form a hard binary, as discussed above. A break ap- 
pears in the profiles at t ss th where the outer, p ~ r~ 2 
profile turns over to a shallower inner dependence; the in- 
ner profile is well approximated as a power-law as well, 
with slope dlogp/dlogr a —1 that gradually decreases 
with time. In projection, this weak power-law cusp pro- 
duces a core-like profile with continuously varying slope. 
Hence this galaxy would be classified as a "core galaxy" 
for t ;> th (Lauer et al., 1995); we note that core galaxies 



also s how weak po wer-law cusps on deprojection (Merritt 
& Fri iman, 1996). We defined the "break radius" Rb as 
the radius where the second derivative of S(-ff) on a log- 
log plot reaches a minimum; this definition is consistent 
with the more common one based on fitting o f a paramct 



ric form to the surface brightness profile (e.g. Lauer et al 



(1995()). Wedefmed 

rb in the same way, as the break radius 

corresponding to the space density profile p(r). Values of 
Rb and rb at several different times are given in Table §. 

Our simulations are the first to demonstrate that weak, 
power-law cusps - corresponding to what are commonly 
called "core" or "cuspy-core" galaxies - can be generated 
by the merger of galaxies with steep cusps, or "power-law" 
galaxies. Since core galaxies are systematically brighter 
than power-law galaxies, it is natural to suppose that weak 
cusps have their origin in mergers. Wc will explore this hy- 
pothesis in more detail below, after discussing the further 
evolution of the density profiles that takes place as the 



BH binary slowly decays. Here we discuss one problem 
with the hypothesis, and a possible resolution. Even some 
moderately bright elliptical galaxies (My w —22) exhibit 
steep cusps, even though these galaxies have certainly ex- 
perienced mergers in the past. How did these galaxies 
avoid the rapid cusp destruction that takes place in our 
simulations? 

A possible answer is suggested by Figures || and (?]. Im- 
mediately after the merger, at f ~ th, the density profile 
is briefly almost homologous with the initial profile, with 
a steep, p ~ r~ 2 cusp. If some mechanism could induce a 
rapid coalescence of the BHs at this time, before the BH bi- 
nary was able to exchange energy with the stars, the steep 
cusp might avoid disruption. We tested this idea using our 
runs in which the two BHs were artificially combined into 
one at t = to + 0.4 = 11.0. The test was successful; the 
density profile after coalescence of the binary (shown as 
the heavy line in the bottom panels of Figure 0) is indeed 
very close to the initial profile and remains so indefinitely. 
We discuss below (jj^j) whether any mechanism might exist 
for inducing such a rapid coalescence, and why it should 
preferentially be active in low-luminosity galaxies. 

4. EVOLUTION OF THE BLACK-HOLE BINARY AND ITS 
EFFECT ON THE STRUCTURE OF THE NUCLEUS 

4.1. Physical Processes 

The evolution of a massive BH binary in a galactic nu- 
cleus has been discussed by a number of authors. We be- 
gin by summarizing that work here and listing the physical 
processes that govern the evolution of the binary and its 
effect on the surrounding stars. 

1. Hardening of the binary. Stars that pass within 
a distance ~ a of the binary, with a the binary's semi- 
major axis, experience a gravitational slingshot and are 
ejected with veloc ities v e j ss Vu n = \J GM\i /a (e.g. Hills 
fc Fullcrton (1980| )); Vbi n is the relative velocity of the two 
BHs if their orbit is circular and Mi 2 = Mi + M2 is the 
total mass of the binary. In a fixed stellar background, 
this leads to hardening at a rate 



d 
~dl 



Gp 



H 



and the rate of energy extraction from the binary is 
dE b G 2 Ml 2 p . 



dt 



-H 



(14) 



(15) 



with H a dimensionless hardening rate. Here p and cr* 
are the mass density and ID velocity dispersion of the 
stars. Equations ( |l4|) and ([Hf ) are derived from a model 
in which the stellar density is assumed uniform and the 
gravitational field from the stars is ignored; gravitational 
focusing by the BH binary is incorporated but not the 
influence of the stellar potential on stellar orbits. The di- 
mensionless hardening rate H is a function of the hardness 
of the binary, measured for instance by Vun/ cr *i as weu as 
the binary's mass ratio M1/M2 and eccentricity. For an 
equal-mass, circular-orbit binary, H varies from ~ 15 for 
an infinit ely hard binary to ~ 2.0 for Vbin I = 1 (Q 
lan, 1996| ). 



Stellar encounters also modify the binary's orbital ec- 
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centricity e. The eccentricity growth rate, 
K ''< 



d\u{i/ a y 



(16) 



is negligible for Vbi n /a ~ 1 and increases to a maximum of 
~ 0.2 for an equal-mass binary with e f=a 0.7 and Vbi n /o~ ;> 
20 flMikkola fc Valtonen, 19921 ; |Quinlan, 1996| ). Changes 
in eccentricity are potentially important because the rate 
of orbital energy loss due to gravitational radiation grows 
steeply for e — * 1 (cf. equation [2^), hence an eccentric 
binary will coalesce sooner than a circular one with the 
same a. 

2. Mass ejection. The binary ejects mass at a rate 



J 



1 



dM P 



M12 dln(l/o) 



(17) 



where J s» 1 is nearly independent of (M1/M2, a) for 
q ;C a/,, and drops with dec r easing hardnes s of the binary 
flMikkola fc Valtonen, 1992| ; |Quinlan, 199^ ). Ignoring the 
variation of J with Vbi n , one can integrate equation (17) 
to obtain 



M ej w JMi 2 ln(q ei /q) 



(18) 



where it has been assumed that mass ejection begins when 
a = a e j] we expect a e j « c%. Thus the binary ejects of 
order its own mass in shrinking from a — a e j to a = a e j/2. 
If the binary's mass is not negligible compared with the 
mass of the pre-existing nucleus, the stellar density near 
the binary will drop as the decay proceeds, causing the 
hardening rate to also drop (equation 14). 

3. Br ownian motion. The binary exhibits Brownian mo- 
tion due to momentum imparted by encounters with stars. 
A single particle of mass M in statistical equilibrium with 
a Maxwellian field of light scatterers with masses m* <C M 
will exhibit an average speed determined by cquipartition 
of energy, 

M 12 (v 2 ) = 3m*al (19) 
and its radius of wandering r w will be given by 



Gp 



(20) 



where p is a mean density averaged over the wandering 
region. These relations ignore any reaction of the back- 
ground particles to the motion of the massive object. Cor- 
rections also apply if the massive object is a binary, which 
receives larger kicks than a point mass from ejected stars. 
The speedup is at most a factor of ~ 2 for a v ery hard 



binary in a nucleus with a steep density profile ( Merritt, 
2001|}J Brownian motion is potentially important because 
it allows the binary to interact with a larger pool of stars 
than if it were fixed at the center of the potential, thus 
prolonging its decay. However the amplitude of the wan- 
dering in an TV-body simulation is likely to be much larger 
than in a real galaxy due to the unphysically small value 
of M12/TO* in the simulations. Brownian motion may also 
help to scatter stars into the binary's sphere of influence 
by introducing a complex time dependence into the grav- 
itational potential felt by the stars. 

4. Loss-cone refilling. Eventually the binary will eject 
most or all of the stars which can come within a distance 



~ q of it. If the binary wanders over a distance r w > a, this 
will happen when it has ejected all stars whose pericenters 
lie within a distance ~ r w from the galaxy center. Once 
the density of stars in the vicinity of the binary drops 
to zero, the binary's decay will stall, unless some process 
can refill the "loss cone." One possibility is infall of a 
third BH, gas cloud, dwarf galaxy or other massive object 
that can perturb the stellar orbits. In the absence of such 
dramatic events, ordinary two-body relaxation will scatter 
stars into the binary's sphere of influence. The associated 
feeding rate is 



dAL, 



M(r max ) 



dt 



(21) 



where M(r) is the stellar mass within r and t r is the star- 



star relaxation time; 



is the radius at which the rate 



of s cattering into the loss cone peaks, typically of order 
r gr ( [Shapiro, 1985 ). 

5. Gravitational radiation. If the decay of the binary 
continues sufficiently far, emission of gravitational radia- 
tion will eventually become the dominant source of energy 
loss. The gravitational radiation time scale t gr is 



t gr — 



gr 64G3 M 3 a 



F(e) 



(22) 



where the factor F(e) contains the eccentricity depen- 
dence: 

F (z) = ; . 73,, ' wr- A (23) 



1 



73.2 
24 c 



37 p 4 
96 c 



(Peters, 1964). The dependence of t gr on e is weak for 



small e; F(0) = 1 and F(0.5) » 0.205. The decay rate from 
gravity wave emission matches that from stellar ejection 
when 



64 G 2 M\ 2 a 



5FH 



c 5 p 



(24) 



The right hand side of this expression is difficult to eval- 
uate ab initio since p will be strongly affected by stellar 
ejection during the binary's decay. However semi-analytic 
models for the combined evolution of the binary and 
the nucleu s (tVIerritt, 2000 ) suggest that a crit ~ 10 -2 'ah- 
Equation (|18|) then implies that the binary must eject 
roughly four times its mass in stars in order to achieve 
gravitational radiation coalescence. We do not include an 
energy sink term corresponding to gravitational radiation 
in our simulations but use equation (|24|) to estimate when 
coalescence would occur. 

A crucial question when interpreting TV-body simula- 
tions is the dependence of the results on N. Real galaxies 
have nuclei with TV « 10 7 , much greater than the particle 
numbers amenable to computer simulation. Fortunately, 
the two processes that most directly affect the evolution 
of a BH binary in a galactic nucleus - hardening and mass 
ejection - have rates that depend only on the local density 
of stars, not on their masses. The rates at which these two 
processes occur in our simulations should therefore reflect 
their rates in galaxies whose mass distributions are similar 
to those in our models. 

However both the Brownian motion of the binary and 
the refilling of the binary's loss cone by two-body encoun- 
ters are TV-dependent processes, and their importance in 
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Fig. 8. — Evolution of the BH binary and response of the stellar environment. Red curves: run A8 (Af./m* = 164); green curves: run A4 
(M,/m t = 328); blue curves: run A2 (M./m* = 655). (a) Separation between the BHs; (b) inverse semimajor axis a -1 = —2Ef,/GM^; (c) 
orbital eccentricity. When calculating a and e before t^, we corrected the BH masses by adding the mass in stars bound to each of them, (d) 
Average density of stars within r < 0.04 from either of the BHs; (e) velocity dispersion in a region 0.01 < r < 0.04 from either of the BHs; 
(f) total mass of stars ejected by the binary, in units of the binary mass. 



our simulations is expected to be much greater than in real 
galaxies. The binary's wandering radius scales as 



-St 



A,r- 1/2 



while the rate of scattering of stars into the loss cone varies 
as 

dM scat 



dt 



oc t, 1 cx cx N 1 



(26) 



where the ./V-dependence of the Coulomb logarithm has 
been ignored. These "second order" effects do not directly 
influence the binary's evolution but they do determine how 
large a supply of stars is available to the binary and hence 
how long its orbit can continue to decay. For instance, a 
wandering binary can interact with a larger pool of stars 
than a binary that is stationary. 

In a real galactic nucleus where TV is very large, Brow- 
nian motion and two-body relaxation would be expected 
to be almost negligible. (Exceptions might occur in nu- 
clei where the gravitational potential is very lumpy, due 



to giant molecular clouds, star clusters, additional massive 
BHs etc.) An obvious inference , drawn by s everal authors 
(IValtonen, 1996| ; pVlerritt, 200"o| |Zier, 20001 ; jGould k RE 
2000), is that a massive BH binary should rapidly eject 



those stars whose orbits bring them within its sphere of 
influence, after which the binary separation should cease 
to change. These arguments are not air-tight however be- 
cause of the complicated way in which the various physi- 
cal processes interact. For instance, mass ejection lowers 
the density of stars, causing the hardening rate to drop 
(equation |l4| ), but the reduction in density leads to an in- 
crease in the wandering radius (equation Gu) allowing the 
binary to move out of the low-density region into a region 
of higher density where it can continue interacting with 
stars. 

4.2. Hardening rate and mass ejection 

Figure || summarizes the evolution of the BH binary in 
our simulations. At any given time, the binary can be 
described by its semimajor axis a, its eccentricity e, the 
direction of its orbital angular momentum vector n, and 
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Table 3 

Measured TV-Body Parameters 



Quantity 


Time 


A8 


A4 


A2 




11 


0.038 


0.036 


0.035 


n 


15 


0.049 


0.048 


0.035 


rt 


19 


0.087 


0.086 


0.056 


Rb 


11 


0.036 


0.033 


0.026 


Rb 


15 


0.053 


0.068 


0.030 


Rb 


19 


0.080 


0.072 


0.063 


da- 1 /dt 


[11, 14] 


970 


920 


750 


da~ 1 /dt 


[16, 19] 


660 


540 


690 


H 


13 


9.3 


7.9 


5.8 


H 


17 


8.1 


6.3 


6.8 


K 


> 11.6 


0.70 


0.13 


0.13 


a e j 


> 12 


0.00098 


0.00085 


0.00074 


a e j 


> 13 


0.00085 


0.00074 


0.00070 


j 


> 14 


0.00090 


0.00071 


0.00070 


J oc 


> 12 


0.45 


0.45 


0.46 


J oo 


> 13 


0.49 


0.49 


0.48 


•I OG 


> 14 


0.48 


0.51 


0.48 


012 


13 


0.087 


0.064 


0.034 




16 


0.085 


0.061 


0.032 


012 


19 


0.081 


0.058 


0.030 


Tid 


[11,19] 


0.028 


0.011 


0.0084 



the position and velocity of its center of mass. The bi- 
nary's hardening rate da~ 1 /dt (Figure |8)b) is nearly con- 
stant with time following the "knee" at t « th when the 
binary first becomes hard. Minute fluctuations in 1/a re- 
flect perturbations in the binary's binding energy due to 
stars tightly bound to one of the BHs; sudden jumps indi- 
cate times when the binary ejects a single star. At low TV, 
the discreteness of individual ejection events induces sta- 
tistical fluctuations in the value of 1/a that in our opinion 
are responsible for most of the ~ 20% differences in 1/a 
between runs with different TV. Table || gives values of 
da -1 /dt obtained by fitting straight lines to 1/a over in- 
tervals 1 < t — th < 4 and 5 < t — th < 8. There is an 
apparent, though slight, decreasing trend of the harden- 
ing rate with TV in the former interval, while the latter 
interval shows no identifiable trend. We also give in Table 
|^ the dimensionless hardening rate H computed from the 
measured values of da -1 /dt using equation (fl4|); the stel- 
lar density and velocity dispersion in that expression were 
evaluated by averaging inside a sphere of radius r = 0.04 
around the binary (for evaluating the velocity dispersion, 
we excluded the center r < 0.01 where stars are strongly 
perturbed by the binary) . We find that H ranges from ~ 6 
to ~ 9, consistent with the moderately-hard binary results 
of three-body scattering experiments summarized above. 
The lack of a noticeable TV-dependence in the hardening 
rate is consistent with the fact that the central stellar den- 
sity and velocity dispersion also do not vary substantially 
with TV (Figure |d,e). 

The orbital eccentricity of the BH binary evolves with 
time as well (Figure |]c:). The eccentricity was evalu- 
ated from the binding energy and the angular momen- 
tum of the binary; before th, we corrected the BH masses 
by adding the mass in stars bound to each of them, 
M, + M*(r < a/2). The initial orbital eccentricity of the 
galaxies is ~ 0.75 and this eccentricity is reflected in the 
relative orbit of the two BHs when they first fall to the 



center. However the orbit rapidly circularizes due to the 
strong density gradients in the cusp, and by t — th the ec- 
centricity is essentially zero. From that point on, e grows 
at an approximately constant rate, albeit with substantial 
fluctuations. The fluctuations do not seem to be strongly 
correlated with TV. Final values of e range from ~ 0.15 to 
^0.3. Table || gives average values of the dimensionless ec- 
centricity growth rate K obtained by fitting the integrated 
form of equation (|l6|), assuming constant K, to data for 
t > 11.6. In runs A2 and A4, the growth of eccentricity is 
well approximated by the relation 



If In 



&ecc \ 
a J 



(27) 



with K — 0.1 3 and a err — 0. 001. (The three-body scatter- 



ing results of |Quinlan (1996 ) predict a substantially lower 
rate, K sw 0.04 for e = 0.3 and Hin/c* = 10-0, although 
they allow for If = 0.13 when e J> 0.4 and Vun/a* ^ 30.0.) 
If the eccentricity continued growing at the rate predicted 
by equation (p7|), it would reach e = 0.5 for aw 2 x 10~ 5 , 
which is near the semimajor axis when the hardening rate 
due to gravity wave emission becomes larger than that due 
to stellar ejection (§||). Equation ( p2] ) would then imply a 
gravitational radiation time scale ~ 5 times shorter than 
if the binary were circular. 

Figure ||f shows the mass ejected by the binary M e j as 
a function of time. We monitored the mass ejection by 
counting stars with positive energies, or cquivalently, with 
galactic escape velocities. This conservative criterion un- 
derestimates the number of stars ejected by a moderately 
hard binary (Vbi n /<J* < 4) because some of these stars are 
not energetic enough to escape the galaxy. Ejected mass 
at the end of each run is close to Mi 2, the total mass of 
the binary. We fit M e j(t) to the integrated form of the 
mass ejection law, equation (p"8|), to derive J and a e j- The 
fit was satisfactory for t ^ th + 2; the fitted values of J 
and a e j are given in Table H. We find J » 0.5 which is 
consistent with the scattering experiments cited above. At 
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earlier times the mass ejection law with constant J under- 
estimates M e j/M 12 . As was the case for the hardening 
rate, the ejected mass does not show a strong dependence 
on particle number in these simulations. 

Readers should not compare the ejecte d mass in our sim 



illatio ns wit h that in the simulations of Quinlan & Hern- 
quist (1997); these authors measured greater M ej - for a 



smaller increase in 1/a, but they also used a different, 
more liberal definition for M e j . The primary focus of our 
study is to relate the dynamics of the BH binary's harden- 
ing to the observable responses of the stellar nucleus, such 
as the decrease in slope of the central cusp, and the precise 
definition of M e j is not important for any of the physical 
conclusions drawn here. 

Although the supply of stars to the binary BH remains 
high throughout our simulations, there is nevertheless a 
steady drop in the stellar density as stars are ejected from 
the core. This can be seen, indirectly, in the slight curva- 
ture of the hardening rate plot (Figure |^b), and directly 
in the change in density within a sphere of radius 0.04 
centered on the binary (Figure ||d) . More detailed infor- 
mation about the change in the stellar mass distribution 
is given in Figure fj]. Most notable is the drop in central 
density between the first profile, at t = th, and the second 
that is offset in time by one A-body unit. This is the same 
rapid drop in density that was discussed above (§2), asso- 
ciated wrth_JonnaU£n_cd^_thc_^^ 



mass sjected during this short time interval is somewhat 
greater with smaller A, which is also evident in Figure 
||d. We attribute this mild A-dependence to statistical 
fluctuations, and also to a spurious effect associated with 
the wandering of the BH binary: since density profiles in 
Figure [t] are centered on the binary, densities at radii less 
than the wandering radius may be artificially lowered. 

There is no suggestion that a "hole" is forming around 
the BH binary; apparently the supply of stars is great 
enough that the binary can eject of order its own mass 
without driving the central density to a very small value. 
The central density profile, which is slightly steeper than 
p ~ r _1 at t = th, becomes slightly flatter than by 
the end of the integrations (Figure [71) and would presum- 
ably become ever flatter if the integrations were continued 
to longer times. The inner density profile remains well 
described as a power law at all times. 

4.3. Brownian motion 

Figure || shows the trajectory of the two BHs at high 
spatial resolution in our simulations. The Brownian mo- 
tion is apparent as a sudden change in the character of the 
BH orbits at t « th- Prior to this time the trajectories are 
smooth and symmetric, reflecting the dynamical-friction 
induced coalescence of the two cusps. However starting at 
t « th the motion becomes more chaotic, resembling a ran- 
dom walk. Figure ^| shows clearly that the amplitude of the 
random motion is a decreasing function of A, as expected 
from equipartition arguments (equation |l^) . We quanti- 
fied the iV-dependence by computing u\i = \J (v 2 ) /3 for 
the binary's center of mass. This quantity exhibits almost 
no evolution with time (Table 0) ; Figure fol shows aver- 

1/2 

ages for t > th- The equipartition relation o~yi °c M 12 
is approximately satisfied. Of interest is the amplitude of 
oi2 which is expected to be slightly larger for a binary 



BH than for a single BH ( Merritt, 2001 ). We were unable 
to check this prediction by a direct comparison between 
our two- and single-BH runs since the latter were not ex- 
tended long enough that an accurate characterization of 
the BH's Brownian motion could be obtained. Instead, 
we compared CTi2 with the velocity dispersion predicted 
by equation (|19|), <7i2 = (m^/M^) 1 ' 2 ^. This requires a 
choice about how to evaluate c* , which depends weakly on 
radius. In Figure [lO] we plot a range of predictions for o~i2 
based on measured values for a* within r = 0.2 (<t* ~ 0.6) 
and r — 0.01 (tr* w 1.0). The rms velocity of the binary 
appears to be slightly greater than expected for a point 
mass, as predicted by Merritt (2001 ). 

Brownian motion of a massive binary was descri bed in 



a few earlier studies. Quinlan fc Hcrnquist (1997 ), in a 



series of ./V-body simulations, noticed a wandering of their 
BH binary with an amplitude 5 — 10 times greater than 
expected on the basis of an equation like (p0|); they at- 
tribut ed the discrepa ncy to inelastic scattering of ejected 



stars. Makino (1997) carried out A-body sim ulations sim- 
ilar to those of Quinlan & Hcrnquist (1997) but using a 
more conservative, direct-summation code and no mass 
spectrum for the field stars. Makino's Figure 7 shows a 
wandering amplitude that scales as ~ TV" 1 ' 2 , and the rms 
velocity of the binary appears to be comparable to that ex- 
pected for a point mass. Makino's results seem consi stent 



with those obtained here and with the predictions of Mer- 



ritt (2001): a massive binary at the center of a dense cusp 



should exhibit only slightly greater Brownian motion than 
a po int particle of the same tota l mass. It is not clear 
why Quinlan & Hcrnquist (1997) found a much greater 
amplitude for the Brownian motion in their simulations; 
some possible reasons for the discrepancy are discussed in 



12 



0.1 



0.05 



0.02 



0.01 



~i 1 r 




200 



500 1000 2000 



M 12 /m, 



Fig. 10. — ID rms center-of-mass velocity of the BH binary's 
Brownian motion, in the merger plane (filled circles) and perpen- 
dicular to the merger plane (empty circles). Equipartition of energy 
implies o\i = (Mi^/m,) -1 / 2 ^ which depends on the region over 
which we average cr*. Solid lines show the band of values consis- 
tent with equipartion between cr* = 0.6 (r < 0.2) and <r* = 1.0 
(r < 0.01), the binary's radius of influence. 
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Fig. 9. — Trajectories of the BHs before formation of a hard binary (blue curves) and after (black curves). The character of the motion 
changes suddenly to a random walk, characteristic of Brownian motion, after the binary forms. The amplitude of the Brownian motion 
decreases with increasing BH mass (as indicated on the left). 
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Fig. 11. — Angular inclination of the BH binary's axis of ro- 
tation. Red dots: run A8 (M./m, = 164); green dots: run A4 
(M./m* = 328); blue dots: run A2 {M./m, = 655). Binary's an- 
gular momentum unit vector in polar coordinates h{9,4>) is shown 
as a dot at position (9 cos <j>, 9 sin</)). Center of the plot corresponds 
to rotation in the merger plane. Inclination of the binary undergoes 
a random walk with an amplitude that decreases with increasing 
M./m,. 



Mcrritt (gOOlp . 

The orientation of the binary is also affected by encoun- 
ters. Figure [ll| shows the direction of the binary's angular 
momentum vector n(9, <p) where each dot represents the 
angular tilt of n from the merger axis. Points shown cor- 
respond to t > th] before the binary is hard, its orientation 
changes at much higher rate, albeit with similar amplitude, 
due to transient bulk perturbations occurring during the 
galactic merger. After th, the bulk torques are negligible, 
but the orientation still changes due to torques imparted 
by elastic and inelastic encounters with stars. It is evi- 
dent in Figure [ll] that the net effect of stellar ejections is 
a random walk in the tilt-space (9, <p) and that the rate 
of orientation-changing decreases both with hardness and 
with N. We observed a maximum tilt of 9 « 12° degrees 
from the merger axis. 

4.4. N -dependence of the evolution 

As discussed above, we do not observe an appreciable 
dependence of the binary hardening rate on N in our sim- 
ulations. This is reasonable since the expected hardening 
rate (equation |l4| ) depends only the mean density p and 
velocity dispersion a, of the field star distribution, not on 
the masses of field stars. 

Some ear lier studies noticed a more ap preciable N- de- 
pendence. Quinlan fc Hernquist (1997 ) found that the 
decay rate dropped with TV until N ~ 10 5 , then seemed 
to level off at N = 2 x 10 5 . |Makino (1997| ) described the 
TV-dependence of the decay rate as weaker than A -1 but 
gave no further details. In both of these studies, the A- 
dependence of the hardening rate was attributed to the 



Brownian motion, since larger N implies a reduced ampli- 
tude of wandering and hence a smaller pool of stars that 
can interact closely with the binary. Quinlan & Hernquist 
showed that the binary's hardening rate dropped rapidly 
to zero if the binary was artificiall y fixed in space. This 
result is also implicit in the study of ^icr (2000| ) who calcu- 
lated the rate of depletion of stars around a massive binary 
that was fixed in space. 

Why do we fail to see any clear A-dependence of the 
hardening rate in our simulations - given that the Brown- 
ian motion does vary in the expected way with N (Figure 
10)? The main reason, we believe, is our very different 
initial conditions, which guarantee a larger supply of stars 
to the binary t han in earlier st udies. We noted above that 
the m odels of Makino (1997) and Quinlan & Hernquist 
(1997| ) had much lower central densities than ours at the 



time of formation of the hard binary. The supply of stars 
was correspondingly smaller, implying a more rapid de- 
pletion of the "loss cone;" once this occurs, the supply of 
stars is essentially cut off and any further decay depends 
on A-dependent processes such as loss-cone refilling and 
Brownian motion. The origin of the muc h lower initial den- 
sity in the Quinlan & Hernquist (1997) simulations may 
be seen in their Figure lc, the evolution of the Lagrangian 
radii for a run where the BHs have a combined mass that 
is 1% of the galaxy mass, comparable to the value in our 
simulations. The stellar mass within a sphere of radius 
0.01 drops from ~ 1% initially to ~ 0.05% by the time of 
formation of the hard binary. In our simulations, the drop 
is from ~ 1% to only ~ 0.3% over a comparable interval 
of time (Figure |). The reason for the much greater den- 
sity drop in the Quinlan & Hernquist simulations was their 
choice of initial conditions: the BHs were initially placed 
far outside the central cusp and fell in. Makino's (1997) 
models had large cores from the start. 

The "supply" of stars in our models can be defined, very 
approximately, as the number of stars (at t = th, say) with 
pericenters less than some critical value p cr it ~ a. This 
definition ignores changes in the stellar density profile that 
result from the changing potential, as well as any back- 
reaction of the binary's motion on the stellar distribution. 
We also ignore any dependence of p cr it on stellar velocity, 
even though low-velocity stars have a larger cross section 
for interaction with the binary than high-velocity stars. 
Figure [l^ shows M cr it , the mass in stars with pericenters 
below p C rit, as a function of p cr it- This was computed by 
counting orbits penetrating, or entirely contained within, 
the sphere of radius equal to p C rit in an N = 2 18 particle 
realization of our initial model (equation ||) of pre-merger 
galaxies. For p cr it <J 0.1, the function is approximately a 
power law 

M / \ 0.84 

Merit -, n I Pcrit \ /'oo\ 

(28) 



M 



1.8 x 



where M is the mass of the galaxy and ro is the half- 
mass radius. The figure also shows l/a cr it, the inverse 
semimajor axis attainable for an ejected mass of M cr it- 
This was computed from the relation 



t/JM 12 



(29) 



using the fitted values of J and a e j (Table ||). Taken to- 
gether, equations (E8h and (E9h relate the size of the region 
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Fig. 12. — Relation between the ejected mass and the shrinking 
of the BH binary. Upper panel is the mass in stars with pericenter 
distances smaller than a given value p cr it as function of p cr a (solid 
curve). For comparison, we show the mass in stars inside radius 
r = Pcrit (dashed curve). M cr ; t is the mass accessible to the BH 
binary provided that it can visit every point within a distance p cr a 
from the center of the galaxy. Lower panel is the maximum inverse 
semimajor axis l/a cr i t at which the binary will stall if it can interact 
only with stars having pericenter distances less than p cr it (solid 
curve) . 



accessible to the binary (p cr it) to the minimum semimajor 
axis that the binary can reach by ejecting all stars visiting 
the region (a crit ). 

In order for the binary to have access to sufficient mass 
allowing it to decay to a -1 w 10 4 , roughly the final value 
in our simulations, Figure [l^ suggests that it must eject 
all stars with pericenters lying within ~ 0.003. This is 
comfortably smaller than the wandering radius even in 
our largest-A^ simulation [r w ~ 0.01; Table ||). We believe 
that this explains why none of our binaries have managed 
to "deplete the loss cone" and stall - their supply of stars 
was never depleted, at least over the interval of time repre- 
sented by the simulations. Combining Figures [lO] and O, 
we estimate that N would have to be increased to ~ 5 x 10 s 
per galaxy in order for the Brownian motion to be small 
enough that the supply of stars would be exhausted by the 
time that 1/a — 10 4 . 

Given that the "loss cone" is never depleted in our sim- 
ulations, how great is the effect of two-body relaxation on 
the decay rate of the binary? We estimated, using stan- 



dard expressions like equation ( pT| ) , that scattering of stars 
into the loss cone probably does make a significant contri- 
bution to the decay rate of the binary in our lowest- iV run. 
For the larger N runs it is probably of negligible impor- 
tance. This conclusion would presumably change if the 
loss cone ever became fully depopulated. 

What would happen in this case - if the supply of stars 
to the binary were depleted, either by continuing the inte- 
grations to much later times, or by using a larger N and 
thereby reducing the amplitude of the binary's Brownian 
motion? Our guess is that the decay would in fact stall, 
producing a BH binary whose separation remained nearly 
constant for extended periods of time. The same would 
presumably also result from initial conditions with much 
lower central density, for instance, the merger of two gi- 
ant elliptical galaxies. We return to the question of the 
persistence of BH binaries in §|[ 

We caution that our simple picture, of a BH binary wan- 
dering against a fixed distribution of stars, is not com- 
pletely correct. In fact the binary "carries" the cusp with 
it to a certain extent. This is clear from the density pro- 
files plot, Figure [?]: the cusp continues as a power law to 
radii much smaller than the wandering radius. It follows 
that the wandering of the binary is a complex process in- 
volving time-dependent changes in the stellar distribution, 
and these changes probably affect to some extent the sup- 
ply of stars to the binary in a way not reproduced in our 
simple analysis. Brownian motion of a massive object in 
a density cusp would be a fruitful topic for further study. 

5. KINEMATICS 

Just as the BH binary affects the stellar density pro- 
file at distances as large as the break radius, so we expect 
the presence of BHs to shape the remnant's kinematical 
properties well beyond the binary's gravitational radius of 
influence. Here we present simulated "observations" of our 
largest-TV merger remnants, from runs A2 and B2, in zero 
dimensions (circular apertures), one dimension (slit) and 
two dimensions (integral field) . The purpose is to generate 
predictions that can be tested with the current generation 
of high-spatial- resolution spectrographs such as STIS, OA- 
SIS and SAURON. Following this, we present several re- 
vealing views of the remnant in ways that are not directly 
accessible to astronomical observation. 

We "observed" our galaxies in two steps (Appendix C|). 
Starting from the stellar velocities projected into an edge- 
on view of the galaxy, we recovered the line-of-sight veloc- 
ity distributions (LOSVDs) non-par ametrically vi a max- 
imum penalized likelihood (MPL) QMerritt, 1997| ). The 



MPL estimate N(V) of an LOSVD is computed on a grid 
in velocity such that it maximizes the log-likelihood of the 
distribution of line-of-sight projected stellar velocities in- 
side an aperture, subject to a penalty function that mea- 
sures the lack of smoothness of N(V). Once N(V) was 
obtained, we expan ded it into its Gauss-Hermite (GH) 
moments de fined by Gerhard (1993[) accor ding to the pre- 
scription of van der Marel & Franx (1993). Of particular 
interest are the four parameters in the GH expansion (Vo, 
a , ha, hi) that quantify, respectively, the mean velocity 
and velocity dispersion of the Gaussian prefactor, and the 
odd and even first-order departures from a Gaussian dis- 
tribution. (Henceforth in this section, the terms "mean ve- 
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Fig. 13. — Line-of-sight velocity distributions (LOSVDs) of the 
merger remnant in the run A2 as a function of aperture diameter 
D. The aperture is circular and centered on the BH binary; its 
diameter varies from D = 0.002 (LOSVD with the shortest peak 
and the broadest wings) to D = 0.2 (LOSVD with the tallest peak 
and the narrowest wings). 



locity" and "velocity dispersion" are used in this restricted 
sense.) These parameters, however, are insensitive to the 
power- law wings expected in LOSVDs in the vicinity of 
a BH ( van der Marel, 1994 ) and hence it is important to 
consider the full LOSVD. 

To increase the resolution inside each circular aperture 
of diameter D, we superposed one hundred snapshots of 
the galaxy that were sampled over 1 ./V-body unit in time, 
corresponding to ~ 10 crossing times at a radius r ~ 0.1 
from the center. Averaging over such a wide time interval 
ensures that stars in the aperture are sampled at random 
orbital phases. In this procedure each dataset was shifted 



in spaice su thai the DII bhiai y lay al Ihe ceiiLei. The cen- 
t ering could cause spurious smoothing on scales 3mallcr 
than the radius of the Brownian motion of the binary (cf. 
§[|) which however amounts to not more than r w 0.0084 
in the run A2 that was used for this purpose. When su- 
perposing datasets, the stellar velocities were left in their 
original frame. Orbits closely bound to the BH binary 
may follow the binary on its random Brownian trajectory, 
thereby incurring a net drift in the velocity. This drift, 
however, scales as (m*/Mi 2 ) 1 / 2 cr* « 0.05 and can be ig- 
nored at this stage. 

We anticipate the discussion in §||.l by quoting typical 
values for the physical scale of our models. Scaling to a 
dwarf elliptical galaxy like M32 gives a factor C « 50 pc 
for converting model dimensions into physical lengths. In 
the case of a giant elliptical galaxy like M87, this factor is 
~ 3 kpc. Note that the radius of gravitational influence 
r gr of the BHs is ~ 0.02, independent of the scaling. 

In Figure [l3| we show a family of LOSVDs from the run 
A2 for a range of diameters D of a circular aperture cen- 
tered on the BH binary. The LOSVD that is most nearly 
Gaussian is seen inside an aperture with diameter about 



D 

Fig. 14. — Fourth Gauss-Hermite moment /14 as a function of 
aperture diameter D. The aperture is centered on the BH binary 
in the run A2 (thick line) and on the single BH in the run B2 (thin 
line). 



twice the BH binary's radius of influence, D ~ 2r gr . The 
LOSVD seen in the smallest diameter aperture D = 0.002 
has the shallowest peak and broadest wings, while the one 
seen in the largest diameter aperture D = 0.2 has the 
steepest peak and almost non-existent wings. These differ- 
ences are reflected in the values of hi shown in Figure ^4[ 
In general, positive /14 indicates that the LOSVD is sharp 
(or "triangular" ) at the top and decays more mildly on 
the sides; negative hi indicates that the LOSVD is broad 
(or "boxy" ) at the top and steep on the sides, reflecting a 
sharp maximum velocity cutoff. 

We observe that hi decreases from hi « 0.12 for D = 
0.003 to hi « —0.024 for D = 0.2 passing t hrough zero for 
D ~ .03. A similar trend was derived by van der Mare] 
(1994|) under the assumptions of isotropy and spherical 



symmetry for a model of M87 with a 5 x 10 9 Mq BH; the 
size of the aperture influences only the overall normaliza- 
tion, and not the velocity dependence, of the LOSVD in 
the large- velocity limit. For comparison, in Figure |lj we 
also plot hi from the run B2 where two BH were combined 
into one at t = t)- t . The absolute amplitude of hi around 
the single BH is smaller than around the binary, which 
is counterintuitive in view of the even wider wings (not 
shown here) that we found for small apertures in the run 
with one BH. The Gauss-Hermite moments, however, are 
sensitive to the velocity profile in the range V <J 2<7o and 
indifferent to the high- velocity behavior. It is therefore 
possible that an LOSVD is both "boxy" at low velocities 
(hi <, 0) and "wingy" in high velocities (which may or 
may not imply a positive hi). 

In Figure O we plot the major axis slit kinematics 
of the merger remnant in the run A2 (BH binary) and 
B2 (single BH). To increase the resolution, we added to- 
gether six views of the major axis, rotated by angles kir/3, 
k = 0, 1, 2, 3, 4, 5, around the minor axis, from the original 
line of sight (k = 0). This yields Vq and /13 that are odd 
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Fig. 15. — Slit kinematics of the run A2 (BH binary; thick line) 
and B2 (single BH; thin line). Slit is positioned along the major axis 
and centered on the BH. The aperture diameter (slit width) varies 
from D = 0.02 at X = ±0.2 to D = 0.004 at the center. Parameters 
Vq, (to, hs and in the Gauss-Hermite expansions of the LOSVDs 
were calculated as described in the text. 



and <7 and h\ that are even under reflection X — > —X. We 
first established the position-dependent maximum aper- 
ture diameters (D = 0.004 at the center and D = 0.02 at 
the ends) and then narrowed the apertures with the re- 
quirement that the number of stars inside each aperture 
remain larger than a fixed amount (4096), ensuring homo- 
geneous statistics across the slit and uncorrelated sampling 
(D < AX) at the center. 

We also computed spatial (unprojected) kinematical 
quantities, in particular, the rotational velocity and 
three diagonal moments of the velocity dispersion ten- 
sor ay , 00 , erg. This was done by averaging inside circu- 
lar wedges of opening angle \8 — 90° | < 30° (major axis) 
and \9 — 90° | > 30° (minor axis). Figure |l^ plots spatial 
kinematical properties of the galaxy immediately follow- 
ing the merger (left column) and of the final galaxy (right 
column) along the major and minor axes. We also pro- 
vide kinematical moments of the merger remnant where 
the two BHs were coalesced into one at t = th (dotted 
line). We will argue below (§||.4) that this model might 
be a good representation kinematically of a "power-law" 
elliptical galaxy. 

All major elements of the slit kinematics that we observe 
in the simulations are also found in the kinematical profiles 
of real galaxies. 

f. Rotation curve. Our models exhibit flat or slightly 



falling rotation curves at the outer radii, \X\ ;> 0.05, with 
Vq/cto ~ 0.5. At these radii, the runs with one and two 
BHs are essentially identical kinematically. At inner radii, 
\X\ <J 0.05, the two rotation curves differ substantially. 
The rotation curve near the single BH rises and appears 
to diverge near the center indicating a nearly-Keplerian 
rotation pattern dominated by the BH's potential. Near 
the binary BH, however, the rotation curve drops, with a 
~ 50% smaller Vq at the knee marking the drop than in 
the case with a single BH. Our choice of the terms "inner" 
and "outer" is not arbitrary; in fact, transition between 
the two regions coincides with the break radius defined 
above (§[|)- We therefore suggest that the same physical 
mechanism responsible for the shallowing of the nuclear 
density cusp also manages to somehow attenuate the cir- 
cumnuclear rotation. 

We propose a mechanism closely related to mass ejection 
by a BH binary (cf. §[|) that leads to precisely this effect. 
When two stars pass near the binary in opposite directions, 
the star on a co-rotating (prograde) orbit is more likely to 
be captured by one of the BHs since it can interact with the 
BH over a larger orbital phase than the counter-rotating 
(retrograde) star. As co-rotating stars are preferentially 
ejected from the nucleus, there will be an increase in the 
relative number of counter-rotating stars, thereby attenu- 
ating the net rotation. In our simulations this effect is too 
small to reverse the direction of rotation in the nucleus 
but it plausibly explains the 50% difference between the 
runs with one and two BHs. This interpretation implies 
a concrete physical prediction: rotation curves in galax- 
ies with shallow central density cusps (or "core" galaxies, 
cf. §|^) should turn over near the break radius and exhibit 
systematically lower rotation to within several dynamical 
radii of the BH than those in galaxies with steep density 
cusps (or "power law" galaxies). 

The rotation curve of our single-BH ( "power-law" ) model 
looks similar to that of M32 as observed with STIS on HST 



(IJoseph et al., 2000[ ). 

2. Velocity dispersions. The central velocity disper- 
sion exhibits a sudden upturn at a distance ~ r gr w 0.02 
from the BHs. The spike is more pronounced in the run 
with one BH {o-Q t , nax — 1.46) than in the run with two 
((To,mai = 1.04). This difference is reduced when the "cor- 
rected" velocity dispersion a = ao(l + y/6h4,) (not shown) 
is used; a is a closer approximation than o~q to the true rms 
velocity ( van der Marel fc Franx, 1993 ). In the run with 



one BH, o~ max — 1.56, while with two BHs o~ rnax = 1.26. 
Nevertheless, in both the uncorrected and corrected ve- 
locity dispersions, we note a systematically lower value 
around the binary BH compared with the single BH out 
to a radius of ~ 0.1. Lower dispersions around the BH 
binary may at first sight strike the reader as unexpected, 
given that the binary injects energy into the system as it 
hardens. But the loss cone is populated largely by radial 
orbits; as the binary captures, ejects, and removes radial 
orbits from the nucleus, the radial dispersion ay drops 

while the tangential dispersion at = \j a % + a e remains 

constant, resulting in a decrease in the average dispersion 



2oJ/3. 



a = - 

We in fact observe tangentially-anisotropic motions out 
to radii r <J 0.03, as shown in Figure [l6| In the bottom 
right panel we show the variation of the anisotropy pa- 
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Fig. 16. — Internal kinematics along the major axis (thick curve) and minor axis (thin curve) immediately after the formation of a hard 
binary (left column) and at the end of simulation (right column). Solid curves: run A2 with a BH binary. Dotted curves: run B2 with single 
BH (major axis only). The major and minor axis moments were computed by direct averaging inside the circular wedges \0 — 90° | < 30° and 

\8 — 90° | > 30°, respectively. First row: rotational velocity v^; thin dashed curve is the circular velocity v c = y/G(M 12 + M(r))/r. Second 
to fourth rows: moments of the rms velocity dispersion (ay, a^, ag). Fifth row: anisotropy parameter (3 = 1 — <t| ja^ where a\ = (<r^ + ct|)/2. 
The tangentially-anisotropic central region (/3 < 0) grows with time, from r rss 0.01 at t = th to r « 0.03 at t = + 7. The simulation with a 
single BH exhibits no appreciable anisotropy. 
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Fig. 17. — Two-dimensional kinematical maps of the merger remnant from run A2. The line of sight is parallel to the plane of the merger; 
the latter projects into the horizontal, or major, axis of the image. Gauss-Hermite parameters (Vo, co, h%, /14) were recovered as described in 
the text (§p|). To increase the resolution inside each aperture, 100 snapshots from 18.1 < t < 19.1 were superposed; this interval amounts to 
> 10 crossing times at radii r < 0.1 from the center. 



rameter (3 with radius in the final models; (3 is defined as 
(3 = 1 - (Jt/crl- In the case of a single BH, < (3 < 0.2 at 
all radii, which we believe is consistent with (3 ~ 0. In con- 
trast, in the binary BH model, —0.4 < j3 < for r < 0.03 
and < (3 < 0.2 for r > 0.03. The anisotropy that we mea- 
sure outside the binary's radius of influence (r o r sa 0.01) is 



consi stent with t hat from the equivalent run of Quinlan & 
Hernc uist (1997 ). These authors detected (3 ~ —1.0 but 



only at radii so close to the center that the calculation of 
(3 ma y depend on the resolution effects and the binary's 



Brownian motion. Zicr (200C), whose binary was fixed in 
space, found even stro nger central anisotropy. We reit- 



that 



erate the prediction of Quinlan & Hcrnquist (1997 
weak density cusps formed by the action of BH binaries 
are tangentially anisotropic. We however disagree with the 
characterization of the anisotropy as "strong;" in fact it is 
mild ((3 <^ 0.5) on scales that can be resolved observation- 
ally. On smaller scales, we warn against the possibility 
of contamination of any inferred anisotropy by stars that 
have recently been ejected or that are interacting strongly 
with the binary. 

3. Third Gauss-Hermite moments. /13 is constant and 
has a sign opposite to that of the velocity parameter Vq ■ 
If /13 and Vq have the same sign, the prograde wing of 
the LOSVD is wider and the retrograde wing is steeper; 
if the signs are opposite, the reverse is true. Figure |l5j 
shows /13 in the runs with one and two BHs; in the for- 



mer \h 3 \ w 0.063, in the latter \h 3 \ « 0.074. We also 
note the average ratios {h 3 ) / (Vq / ctq) ~ —0.12 for one BH 
and (/i3)/(Vo/uo) ~ —0.15 for two. These ratios are in 
excellen t agreement with th e empirical h 3 - a/V rela- 
tion of Bender et al. (1994 ) in a sample of 44 elliptical 
galaxies. All galaxies in the Bender et al. sample show 
opposite signs of h 3 and Vq and fall near the relation 
(/i 3 ) ss — 0.12(Vo/ao). The sudden sign change in I13 very 
near the center is a consequence of averaging ov er aper- 
tures" similar features can be seen in the models of |DchncE 
(19951 ) and |Qian et al. (1995[) 



Interestingly, our si mulations present a direc t counterex- 
ample to the results of Burkcrt & Naab (2001) whose sim- 
ulated mergers always yielded /13 / '(Vq / cr ) > in apparent 
disagreement with the observations; these authors argue 
that the presence of disk-like subcomponents may be nec- 
essary to reproduce the corr ect sign of ^/(Vn/oYi) < 0. 
The bulge initial conditions of Burkert fc Naab (2001 ) had 
p ~ r _1 central cusps and were thus significantly less con- 
centrated than ours. New simulations of mergers with a 
range of density profiles are needed to clarify the depen- 
dence of /13 on initial conditions. 

4. Fourth Gauss-Hermite moments, is very small 
except at the very center. At radii greater than twice 
the BH radius of influence r gr , our data are consistent 
with hi = 0. For r gr <J r <J 2r gr , /14 dips into negative 
values; this is again a consequence of averaging over finite 
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Fig. 18. — Structure of the merger remnant in the meridional plane. The BH binary is at the lower left corner of each panel. Plots are 
averages over multiple snapshots as described for Figure H. Horizontal axis is distance w from the Z axis; vertical axis is distance Z from the 

equatorial plane, p is the density of stars, v^, is the average velocity through the meridional plane, a = \J <tj: /3 + 2<r^ /3 is the average ID 
velocity dispersion, and (3 = 1 — a'f/a'f, is the anisotropy parameter, where ay and ay are, respectively, the radial and the tangential velocity 
dispersions. 



apertures (Dehnen, 1995; Qian et al., 1995). For r <; 



hi becomes positive again with peak values of about 0.027 
(one BH) and 0.083 (two BHs). Although it is tempting to 
ascribe this difference to the differing effects of single and 
dual BHs, the Poisson uncertainties in our determinations 
of /14 are large enough that we are reluctant to infer a 
significant difference between the two models. 

The ne w generation of inte gral field spectrographs like 
SAURON (Peletier et al., 2001) can obtain two-dimensional 
maps of the stellar kinematics from absorption- line spectra 
at resolutions of less than an arcsecond, corresponding to 
<^ 1 pc in nearby galaxies. In Figure we show 2D maps 
of the parameters Vo, ctq, /13 and /14 similar to the maps 
obtainable with SAURON. The line of sight is parallel to 
the plane of the merger, which projects into the horizontal, 
X-axis of the image. To generate these images, we com- 
bined 100 snapshots in the time interval 7.5 < t < 8.5. At 
every point on an 41 x 41 grid we started with an aperture 
of diameter D = 0.01 that was then shrunk to smaller D 
under the condition that the number of stars inside the 
aperture always remain sufficient for the recovery of the 
LOSDVs (a few x 10 3 ). 

The rotation pattern seen in projection is symmetric 
and co-aligned with the merger's initial orbital axis (ver- 
tical axis in the image). The rotation reaches a maximum 
of about I Vb| ~ 0.4 at r « 0.04 on the major axis before 



dropping slightly at larger radii. The velocity dispersion 
Co exhibits no deviation from a circular pattern in spite 
of the noticeably elliptical isophotes of the density (Fig- 
ure |J). The third Gauss-Hermite moment /13 is largest 
along the major axis (/13 ~ ±0.1 at B = 0.05 is typical) 
and gradually decreases to zero toward the minor axis. 
Its map resembles the rotation pattern of Vb, except that 
the sign of ft.3 is opposite from Vq as noted above. The 
fourth Guass-Hermite moment /14 is consistent with zero 
everywhere except for the center; however there is a hint 
of positive /14 along the major axis and negative /14 along 
the minor axis. 

In the run B2 where the BH binary was replaced by 
single BH at t = th (not shown here), the rotation field 
peaks at Vo « ±0.55 much closer to the BH at \X\ < 0.01 
and remains high |Vb| <; 0.5 in a disk-like circumnuclcar 
region out to \X\ ~ 0.025. 

Finally, in Figure |l8| we present 2D maps of p, v^, a 
and /3 in the meridional plane, which is perpendicular to 
the merger plane and contains the axis of rotation of the 
merger remnant. Here a = ^/a%/3 + 2of /3 is the average 
ID velocity dispersion and f3 = l — a^/al is the anisotropy 
parameter. In each panel, the pixels are mapped in coor- 
dinates (zi7, Z), where w is the distance from the axis of 
rotation and Z is the distance above the equatorial plane 
(i.e. the plane of the merger). We found symmetry in all 
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of the above variables under reflection Z — > — Z ', there- 
fore only positive values of Z are plotted. As in Figure 
17, the images were generated from a superposition of 100 
snapshots in the interval 7.5 < t < 8.5, but this time we 
calculated non-parametric kernel estimates for all quanti- 
ties as described in Appendix |b[ 

6. DISCUSSION 
6.1. Scaling 

The scaling of our models to physical units will depend 
on which time step we choose to compare with real galax- 
ies and on which galaxy we choose for comparison. For- 
tunately there are two quantities that remain nearly or 
exactly constant with respect to time in our simulations 
and which are convenient for scaling: the total mass M\i 
of the BH binary, and the velocity dispersion er* of stars in 
the nucleus, outside the region where the stellar motions 
are strongly affected by the BHs. The first quantity is pre- 
cisely constant, while the second varies only slightly with 
time (e.g. Figure ||e) and position (Figure |l6| ). 

Furthermore there is a tight relation between BH mass 
and stellar velocity dispersion in real galaxies, the M, — a 
relation, which allows us to reduce the scaling to a single 
number. For our purposes, the most useful form of the 
M, — a relation is 



The length and time scaling factors are 



M. = 1.30(±0.36) x 1O 8 M 



4.72(±0.36) 



( |Merritt fc Ferrarese, 2001a|) 



200 km s" 

(30) 

Here er c is the projected ve- 
locity dispersion measured in an aperture of radius r e /8 
centered on the BH, with r e the half-light radius of the 
stars. At ground-based resolutions, a c is essentially unaf- 
fected by the presence of the BH and measures the velocity 
dispersion defined by the stellar spheroid. If we equate a c 
with cr* in our definition for ah, the semi-major axis of a 
hard binary (equation |]), we find 



ah 



1.51 pc 



( M. \ 



0.576 



(31) 



Vio 8 aW : 

In our simulations, a*(r e /8) is close to 0.8 in model units 
at all times after formation of the hard BH binary. The 
mass of the BH binary is 0.02 in model units; we identify 
this with M., the mass of the BH in the observed galaxy 
(or the combined mass of the two BHs in the case of a 
binary). If we define scaling factors {M, V, £, T} for our 
models, such that the mass in physical units is M. times 
the mass in model units and similarly for velocity, length 
and time, then ( pp| ) implies 



0.02.M = 1.30 x 10 8 M, 



or 



M 



2.27 



0.8V 

^ 200 km s" 1 

V 



4.72 



and 



1O 9 M V200 km s^ 1 

M = 50M.. 



(32) 

(33) 
(34) 



390 pc 



/ M. 



V 1O 8 M 



0.58 



T 



1.62 x 10°yr 



/ M. 

r . „ _ _ 



V1O 8 M 



0.37 



(35) 



We note that these scaling factors are independent of quan- 
tities like the break radius r\,. This is appropriate, since 
the empirical M, — a relation is also (apparently) rv in- 
dependent, and in our simulations, (and hence V) are 
hardly affected by the creation of a core. 

We consider two representative examples for the scaling. 
Suppose we identify our simulations at early times (before 
cusp destruction) with a galaxy like M32, a dwarf elliptical 
with a s teep cusp. The BH mass in M32 is M, 3 x 



1O 6 M (Joseph et al., 2000| ), giving 



M = 1.5xl0 8 A%, V = 110kms" 1 , 
C = 51 pc, T = 4.4 x 10 5 yr. (36) 

The elapsed time in our simulations from t ~ th, the time 
of formation of the hard binary, until the final time step 
at t ss 20 then corresponds to At 4.0 x 10 6 yr. The final 
value of a, the semimajor axis of the binary, is ~ 5 x 10 -3 
pc or ~ 0.04a/!, and the final break radius is ~ 3 pc. 

Or we could identify our simulations at late times (after 
cusp destruction) with a galaxy like M87, a bright elliptical 
with a weak cusp. The BH m ass in M87 is M, « 3 x 



1O 9 M ( |Macchetto et al., 1997| ), which gives 



M = 1.50 x 10 n M Q , V = 490 kms _1 , 
C = 2.8 kpc, T = 5.7 x 10 6 yr. (37) 

The corresponding elapsed time is ~ 5.0 x 10 7 yr, the 
final value of a is ~ 0.28 pc (also ~ 0.04et/j), and the 
final break radius is ~ 170 pc. (This value is ~ 3 times 



smalle r than the break radius reported by Lauer et al 
(1992] ), which probably means that M87 has undergone 



more than one major merge r sin ce the era of formation of 
the supermassive BHs; see §6.3.) 

Henceforth we will refer to these as the "M32" and 
"M87" scalings respectively. 

6.2. Cusps and Cores 

Our simulations demonstrate that the merger of two 
galaxies with steep, power-law density cusps (p ~ r~ 2 ) can 
produce a galaxy with a shallow power-law cusp (p ~ r _1 ) 
inside of a break radius rj; the necessary ingredient for the 
transformation is en ergy input from a pair of massive BHs. 
Omitting the BHs (Barnes, 1999), or artificially coalesc- 
ing them immediately after the merger (§g), preserves the 
steep cusp. As discussed above (§||), most of the evolution 
of the central density profile in our simulations takes place 
during a brief period when the two BHs first form a hard 
binary. Subsequent ejection of stars by the BH binary 
produces a gradual flattening of the inner slope (Figure 
|7|). The space density profile p(r) is always well described 
as a power law, p ~ r~ 7 , at small radii. However the pro- 
jected density looks qualitatively different: a "core" 
appears, characterized by a log-log slope that falls to zero 
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r 1 cusp 



at the center (Figure 0) . This is natural, since an 
in space density projects to a logarithmic core in surface 
brightness (e.g. Dchncn (1993 )); only power laws steeper 
than r _1 remain power laws on projection. 

Figure invites comparison with the luminosity profiles 
of real galaxies. These are well known to fall into two 
classes, the "power laws" and the "cores," the latter having 
well- defined break radii. It was initially argued ( Laueq 
et al. J 1995Q that core profiles were qualitatively distinct 
from power-law profiles, but Merritt & Fridman (1996) 
pointed out the effect of projection on a weak power-law 
cusp and used nonparametric deprojection techniques to 
verify a power-law dependen ce of p on r even in th e core 
galaxies, 
et al. 



More recent work (Gebhardt et al., 1996; Rest 



2001) has verified this result in larger samples. Our 



simulations demonstrate that nuclear density profiles with 
a range of power-law slopes can be generated in a natural 
way starting from galaxies with steep cusps, and that the 
resultant nuclei look like classical cores in projection when 
the index of the power law is small. 

Surface brightness data for elliptical galaxies and bulges 
are commonly fit to a parametric model, the "Nuker" law: 



s(i?) = s r i (i+n 



a s(r-/3)/a 



R 
Ro 



(38) 



( Laucr et al., 1995 ; Byun et al., 1996). This functional 



form (and the one due to Ferrarese et al. (1994) on which 
it was based) has a built-in power-law dependence of S 
on R at small radii. Our simulations preserve a power 
law in the space density at small radii but not the surface 
density. We suggest fitting power-law models like (38|) to 
the deprojected density profiles of galaxies, rather than to 
their surface brightness profiles - or equivalcntly, fitting 
the projection of an expression like equation (5sT) to the 
surface brightness data. If "core" galaxies really are char- 
acterized by weak inner power laws in p, as in our models, 
the fit to the data should thereby be improved. 



6.3. The M, — and M, 



M e j Relations 



Part of the motivation for putting core and power- 
law galaxies into distinct categories was the observation 
that core galaxies are systematically more luminous than 
power-law galaxies. Core galaxies have —24 <J My < —20 
while power-law galaxies are mostly fainter than My — 
-20 ([Ferrarese et al., 1994 iGebhardt et al., 1996|; |Fabei 



et al., 1997), although with considerable overlap at inter- 
mediate luminosities, —22 <J My <; —20.5. Recent stud- 
ies (e.g. ICarollo k Stiavelli (19981) ; [Rest et al. (200l| ); 



Ravindranath et al. (2001)) have confirmed this system 



atic difference while weakening the case for a dichotomy; 
the variation of cusp slope with galaxy luminosity is es- 
sentially continuous in the larger samples now available. 

What does our model predict? Bright galaxies should 
have experienced more mergers than faint galaxies and 
suffered more from the scour ing action of binary BHs 



This is the r easoning that led Ebisuzaki, Makino & Oku 



mura (1991) to suggest that the "cores" of giant ellipti- 
cals (which they took to be regions of constant density, 
not weak power-la w cusps) are generated by binary BHs. 



Faber et al. (1997) showed that the "core masses" of bright 



ellipticals scale with galaxy l uminosity in roughly the sam e 
way as in the simulations of Quinlan & Hernquist (1997). 



A point not made by these authors is that a pair of BHs 
will eject of order its combined mass during each merger. 
The total mass ejected depends both on the final BH mass, 
and on the number of stages in the merger hierarchy that 
have occurred since the BHs first formed. If the masses of 
the BHs at some stage in the merger hierarchy is M,/n, 
with M. the final BH mass, the mass ejected by all the pro- 
genitors at this stage is of order n x M./n ss M,, and the 
total mass ejected in the complete set of mergers scales 
both with M t and with the number of mergers. Since 
the latter is bigger for bigger galaxies, we expect to see 
a steeper-than-linear relation between core mass and BH 
mass in elliptical galaxies. 

To sharpen this argument and test it against real galax- 
ies, we need a working definition of "core mass," or more 
precisely, for the mass deficit - the mass ejected by the 
BHs. We define this as the mass needed to bring an ob- 
served density profile to a p ~ r~ 2 dependence near the 
center. Here we are assuming, as above, that r~ 2 cusps 
were universally present before the binary BHs began to 
do their damage and that they would have been preserved 
in the absence of BHs. This assumption (similar to the 
one made in §[j] when justifying our choice of initial con- 
ditions) is reasonable since: (1) the growth of single BHs 
in pre-existing cores produces p ~ r~ 2 cusps; (2) steep 
cusps are preserved during mergers in the absence of en- 
ergy input from supermassive BHs (§||); (3) without BHs, 
pre-existing cores evolve int o something like steep cusp s 
through successive mergers (Makino & Ebisuzaki, 1996); 
(4) faint elliptical galaxies universally have steep cusps. 



Consider then a pair of galaxies with p 



central 



density cusps which merge to form a galaxy with a shal- 



lower cusp, p 



7 < 2, inside of a break radius rb- 



Assume that the density profile of the merger remnant was 
homologous with that of the merging galaxies before the 
BHs began to heat the stars. The mass initially within rb 

was 

2o-ln 



4-7T / dr r -— - 

/„ 27rGr 2 



G 



and after mass ejection, 



4tt 



dr 



l-KGrr 1 ^ 3-7 G 



(39) 



(40) 



We ignore changes in er* (cf. Figure g). The ejected mass 
is therefore 

2(2 - 7) -In (41) 



M, 



7 



G 



The break radius r& that appears in equation ( f4l| ) refers 
to the space density p(r), while published break radii Rb 
are derived from surface brightness profiles £(ii). However 
the d efinitions of both r & and Rb are to an extent arbitrary 
(e.g. Faber et al. (1997)) especially since we ignore in our 



definition of M e j the precise form of the density profile 
near the break radius; furthermore we find that rb ~ Rb 
within the uncertainties in our A-body models (Table [}]) . 
We therefore feel justified in replacing rb by Rb in equation 
(|4l|). We find that M ej - defined in this way is larger by a 
factor (~ 2 than the ejected masses that we measured in our 
simulations by counting stars that completely escape the 
galaxy (Figure ||f). There is a simple explanation for this 
discrepancy: after the ejection of one BH-binary-mass in 
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stars, the central gravitational pull on the remaining stars 
will decrease by roughly a half; these stars will then shift 
to wider orbits thereby increasing the apparent M e j . But 
this factor of ~ 2 is comparable to other uncertainties in 
the definition of M e .j and we will neglect it in what follows. 
' Faber Ct al. (1997 ) list 16 "core" galaxies for which there 
are published measurements of 7, and <r*. Essentially 
none of these galaxies has an accurately-determined BH 
mass, but we can use the M. — a relation (f}0|) combined 
with measured values of 0% to estimate M # . Figure [ll] 
shows the result. The correlation of M e j with M, is rea- 
sonably tight, and the slope of the relation is significantly 
greater than one, as predicted. By contrast, Rb and M, 
are essentially uncorrelated. The brightest core galaxies 
(My S3 -24) have (M ej ) S3 10M. while the faintest core 
galaxies (My S3 -20) have M ej S3 M,. Ejection of ~ 10M. 
in stars in the biggest galaxies seems reasonable since these 
galaxies are believed to have experienced several mergers 



since the quasar epoch (Kauffmann, Chariot & Balogh 
200T| ). 

While this success is encouraging, we point out that 
there are other factors that might contribute to the trend 
of increasing M e j/M, with galaxy mass. The "power-law" 
galaxies have 7 S3 2 and hence M e j S3 according to equa- 
tion (fill). These galaxies would fall far below and to the 



o 
o 
o 



o 
o 



R b (pc) r " 
o 



(a) 



-H+| 1 — I I I I 1 1 1 — I I I I 1 1 1 H Hfat+Hj 



(b) 




10 



10 



M. (M Q ) 



Fig. 19. — Correlation of BH mass with break radius (a) and 
ejected core mass (b) in galaxies classified as "core" galaxies by 

Faber et al. (1997} .BH masses were computed using the M, — a 

relation (equation p0| ) and measured values of the central velocity 
dispersion. M e j is defined as in equation (kit); the inner slope 7 = 
— dlog u/d log r was taken from Gebhardt et al. (1996). Solid line 
in the lower panel is a least-squares fit. Dotted lines show M e j = 
(1,2,4, 8,16, 32) M.. 



left of the "core" galaxies in Figure [Hj and so it is reason- 
able to expect a steeper-than-unit slope for the brighter 
galaxies that are plotted there. We return in the next sec- 
tion to the question of how the power-law galaxies man- 
aged to maintain steep cusps in the presence of mergers. 

6.4. Persistence of Steep Cusps in Faint Galaxies 

The story just outlined can not be complete, since steep 
power-law cu sps persist in elliptical galaxies as bright as 
My —22 (Gebhardt et al., 1996), and are universally 
present in galaxies fainter than My S3 —20. How have the 
cusps in the power-law galaxies managed to avoid destruc- 
tion by merging BHs? We discuss several possibilities. 

1. Power-law galaxies do not contain supermassive BHs. 
While the presence of supermassive compact objects has 
only been reli ably established in a handfu l of galaxies (see 
discussion in Ferrarese & Merritt (2000)), several of the 
best cases are in stellar systems with steep cusps (e.g. the 
bulge of the Milky Way; M32). 

2. Power-law galaxies were not formed via mergers. 
This would contradic t standard models fo r hierarchical 
structure formation ( Lacey fc Cole, 1992 ; Hachnclt & 



Kauffmann, 200C; Menou, Haiman & Narayanan, 2001). 



although it is possible that most power-law galaxies have 
not experienced major mergers since the era when BHs 
gained most of their mass; we return to this idea below. 
3. Cusps in power-law galax ies are regenerated follow- 



Fabcr ct al. (1997j ) suggested that steep cusps 



vng mergers. 

might be produced by star formation from fresh gas sup- 
plied during mergers. While mergers certainly lead to en- 
hanced star formation, we consider this explanation un- 
likely since nuclei formed from infal ling gas do not re- 
sembl e featureless power laws (e.g. Mihos & Hernquist 
|(1994 )). The scale- free nature of the cusps suggests to us 
a gravitational origin. 

4. A mechanism exists for extracting energy from BH 
binaries before they can heat the surrounding stars. This 
suggestion is motivated by the observation that the cen- 
tral density profiles in our simulations remain homologous 
with the initial profile, p ~ r~ 2 , for a short time after the 
merger (Figure 0). If some process were effective at ex- 
tracting the binary's energy at this stage, at a rate higher 
than the hardening rate due to stellar ejection, cusp dis- 
ruption could be avoided, producing a coalesced BH bi- 
nary in a steep cusp. We tested this idea in a limited 
way above (^) by artificially combining the two BHs im- 
mediately after the merger; the profile retained its steep 
power-law character thereafter. 

We propose that explanation (4) is the correct one and 
that the mechanism which extracts energy from the BH bi- 
nary is gas dynamical in origin. The effects of gas on the 
evolution of BH binaries have been discussed by a number 
of authors. Gas may accrete onto the larger of the BHs 
causing orbital contraction at a rate ~ M3 9 / M ( Bcgclman"] 



Blandford & Rees, 1980; Valtoncn, 1996). Interaction of 



the binary with an accretion disk will transfer angular mo- 
mentum f rom the BHs to the gas c au sing the binary orbit 
to decay (Lin fc Papaloizou, 1979a| ,§; [Byer fc Clarke, 1995| 



Artymowicz fc Lubow, 1996[); therate is again of order th e 



gas accretion rate ( [Ivanov, Papaloizou fc Polnarev, 1999 ). 
Discr ete gas clouds, lik e those near the center of the Milky 
Way (Coil & Ho, 1999), might also affect the evolution of 
the binary, particularly when its separation is large. Gas 
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clouds could scatter s tars into the BH binary's loss cone 
(Kim & Morris, 2001), enhance the Brownian motion of 
the binary, or even perturb the individual B H trajectories 
before they form a bound pair ( Bekki, 2000 ). 

A difficulty with explanations that demand gas accre- 
tion is the high accretion rate implied: a mass of order 
Mi 2 must find its way to the nucleus on the time scale 
of the merger, or roughly one crossing time, so that the 
binary avoids disrupting the stellar cusp. Much of the en- 
ergy released from the BH binary would go into heating 
the gas. It follows that the events we are envisioning are 
en ergetically similar to q uasars, a conclusion reached also 
by |Graham ct al. (1990| ), |Gould fc Rix (2000[ ) and others. 
If this picture is correct, power-law galaxies acquired most 
of their mass before and during the quasar epoch, from 
progenitors that were gas-rich, while core galaxies are el- 
lipticals whose most recent major merger occurred after 
the era of formation of the BHs. 

There is circumstantial evidence in support of this pic- 
ture. The division between core and power-law gal axies 
occurs at -22 < My < -20.5 flFaber et al., 1997j ). In 
semi-analytic models for galaxy formation, the predicted 
ratio of gas mass to stellar mass during the last major 
merger is a steep function of galaxy luminosity, dr opping 
from ~ 3 for My = -18 to - 0.3 for My = -21 ( [KamT 
mann| & Hachnclt, 2000). Thus the gas content of the 



progenitors of the power-law galaxies was likely to have 
been high at all previous stages in the merger hierarchy, 
while for core galaxies the most recent mergers were prob- 
ably gas-poor. Furthermore the redshift of the last ma- 
jor merger (defined as a merger with mass ratio less ex- 
treme than 1 : 3) is a strong function of a galaxy's current 
mass. Most galaxies with masses less than ~ 10 10 Mq 
(My « —18) have never experienced a major merger; only 
galaxies with M £ 10 n M Q (My w -21) have typically 



under gone a major merger sin ce a redshift of 1 (Kauff- 
mann 



Chariot fc Balogh, 20010 



Sharpening these arguments will require iV-body simu- 
lations of BH binary evolution including gas. 

6.5. Coalescence Time Scales and Persistence of Binary 

BHs 

We identify two characteristic time scales associated 
with the decay of the BH binary in our simulations. The 
first, which is essentially ./V-independent, is the brief period 
following the merger when the two BHs fall to the center 
and form a hard binary. As discussed above (§||), the BHs 
initially come together in a time that is approximately as 
long as the merger itself; most of the energy transfer from 
the BHs to the stars in the nucleus takes place in just 
At « 0.1 in model units (Figure §), or ~ 10 5 - 10 6 yr. 
The separation between the BHs at the end of this period, 
t w th , is ~ 10~ 3 in model units (Figure [|) corresponding 
to roughly 0.05 pc (M32) or 3 pc (M87). 

The second time scale is associated with the gradual de- 
cay of the BH binary. We especially wish to know how long 
it will take the binary to decay to the point that emission 
of gravitational radiation becomes the dominant energy 
sink. This occurs when |d/a| _1 due to mass ejection first 
equals t gr as defined above (equation [22]). This second 
time scale is potentially A-dependent, since \a/a\ depends 
to some extent on collisional processes which may be much 
larger in our simulations than in real galaxies (§|]). 



For the moment, assume that the decay rate of the BH 
binary in our simulations is characteristic of real galaxies 
with much larger N. The inverse semimajor axis increases 
roughly linearly with time in our runs (Figure |b), which 
suggests that we define a decay rate 



at \ a 



In model units, S w 7.0 x 10 2 ; in physical units, 

-0.95 



(42) 



S w 1.0 x 10 



V 10 8 M. 



yr-^c" 1 . (43) 



We explore the consequences of assuming that S remains 
constant at this value until gravity-wave coalescence takes 
place, then discuss the reasonableness of this assumption. 

Energy loss due to gravitational radiation dominates 
that from stellar interactions when 



tn 



1 

aS 



with t gr given by equation 



64 G A M( 2 
55 c 5 



(44) 



(45) 



and -F(e) (equation g3|) is henceforth set to one. (The mod- 
est rates of growth of e in our simulations, [Q, imply F m 1 
at all but the latest stages of the merger, and a oc F~ - 2 .) 
Combined with equation ([43]), this condition becomes 



a < a c 



0.012 pc 



( M. 



V 1O 8 M 



(46) 



Scaling to real galaxies we find a cr u ~ 8 x 10 -4 pc (M32) 
and ~ 0.2 pc (M87). These values are smaller, by respec- 
tive factors of ~ 0.12 and ~ 0.6, than the final value of a 
in our simulations. The gravitational radiation time scale 
when a — a C riti which is also approximately equal to the 
time elapsed in reaching a cr n from au in this simple model, 



tgr(a cnt ) ps 5.0 x 10 7 yr 



/ M. 



V10 8 M( 







0.16 



(47) 



which is 3.0 x 10 7 yr (M32) and 9.0 x 10 7 yr (M87). These 
are factors of ~ 8.0 (M32) and ~ 1.8 (M87) longer than 
the elapsed time from t — t^ until the final time step in our 
simulations. Thus a straightforward extrapolation of our 
A-body results implies that BH binaries would achieve 
gravitational radiation coalescence in a relatively short 
time, of order 10 8 yr, following a merger. 

Is this a reasonable conclusion? Although the decay 
rate of the BH binaries in our simulations showed no ap- 
preciable A-dependence over the range of A that we tested 
(Figure ||b), other aspects of the evolution were observed 
to depend strongly on A, including the amplitude of the 
Brownian wandering (Figure |9|). We argued above (§|) 
that efficient decay of the binary in our simulations de- 
pended on this wandering, since it allowed the binary to 
interact with a larger pool of stars than if it remained pre- 
cisely fixed at the center. We predicted that the supply of 
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stars would have been exhausted by the end of our sim- 
ulations if the wandering had been reduced by increasing 
N to ~ 5 x 10 5 , still much smaller than in real galaxies. 

In a real galactic nucleus, the random velocity of a su- 
permassive BH is expected to be only 



simulations.) However the basic conclusion is unchanged. 
In galaxies with initially shallower cusps the decay would 
be expected to stall at still larger separations. 

The inefficiency of stellar-dynamical processes at bring- 
ing together superma ssivc BHs has been noted b y a num- 



berof authors (e.g. Polnarcv fc Rccs (1994 ); Valtonen 



1996| ); |Mcrritt (2000| ); |Gould fc Rix (2000| )). We believe 



0.033 km s" 



1/2 



f M. 

M Q J \io 8 m & 



-0.29 that the problem has often been overstated due to the 

(48) use of over-simplified models for describin g the interac- 
tion o f stars with the binary; for instance, Gould & Rix 



of order unity, for a binary BH with M12 = M, (Merritt 
(200l|; Figure |l0|) . ) The corresponding wandering radius 
is 



0.010 pc 



/ M. 



V1O 8 M 



-2.35 



100 pc 



(49) 



the aL - 1 relation has been used to express a, in terms of ( 200C D [ Z nore the J &ct that most of the hardening comes 
M.. (This velocity would be increased by a modes t factor, fr ° m stars on orblts w ' apocenters much greater than 

the semimajor axis a. But the basic argument is sound: a 
fixed binary would have difficulty interacting with several 
times its own mass in stars even if located at the center of 
a steep density cusp, hence it could probably not achieve 
gravitational radiation coalescncc in a time shorter than 
the age of the universe. Unless additional mechanisms ex- 
ist for extracting energy, its decay would be expected to 
stall. 

We discussed above what these "additional mechanisms" 
might be and argued that they would be most effective in 
galaxies whose progenitors were gas-rich. These galaxies, 
which in our opinion are least likely to contain BH bina- 
ries, are also the systems in which detection of BH binaries 
would be easiest, through the measurement of periodically 
varying features in the emi ssion line systems associated 



with r c the core radius inside of which the stellar density is 
taken t o be constant. These relations are based on a simple 
model (Merritt, 2001) in which the stellar distribution is 
assumed to be unaffected by the motion of the BH or BH 
binary. If we accept these figures, the wandering radius of 
a BH binary in a real galaxy is of order or less than the 
semimajor axis of the binary. It follows that Brownian 
motion can generally be neglected when considering the 
inter action of the binary with surrounding stars. 



Thi p conclusion may be too pessimistic. A nearly sta- 
tionary binary would soon "s cour clean" a nearly-spherical 
region of radius ~ 2a — 3a (Zier, 2000), after which the 
force acting on it would be essentially zero and even a very 
small Uui would translate into a large displacement. The 



with one or both of t he BHs (Begelman, Blandford & Rees 



1980; Gaskell, 1995). Such features have tentatively been 



detected in a few active galaxies; the best case is proba- 
bly OJ 287, a blazar with nearly-periodic outbu rsts dating 
back roughly 100 years (Pursimo et al., 200C). However 



the interpretati on in terms of a binary system (Lehto & 



ampli tude of the Brownian motion might therefore slowly Valtonen, 1996| ) is not airtight. A larger number of inter 



increase with time as the binary eats its way through the 
nucleus. Our simulations tell us nothing about the plau- 
sibility of this scenario since our runs never reached the 
point of loss-cone depletion. Larger N, longer runs, or 
initial conditions with lower central densities would be re- 
quire d to test this idea. 



acting systems exhibit emission from two resolved peaks, 
probably the nuclei of galaxi es in the early stages of merg- 



ing (e.g. Carico et al. (199C)); the smallest projected sep- 



But [~suppose that the Brownian motion remains always 
small, with an amplitude less than the separation between 
the BHs. We can conseratively estimate the final hardness 
achieved by the binary by the following argument. Assume 
that the center of mass of the binary remains fixed, and 
that the shrinking binary ejects stars in order of their peri- 
center, from smallest to largest. In this way, the binary 
acts to reduce the central stellar density in the most rapid 
possible way, causing its decay to stall in the minimum 
time. This process would create a hole at the galaxy's 



aration, in ARP 220, is ~ 360 pc ( |5covillc et al", 1998| ). A 
recently-discovered double quasar contains tw o peaks with 
a projected separation of a few kiloparsecs (Junkkarincn 



et al., 2001) 



The supermassive BHs in these interacting systems have 
not yet formed bound pairs. True BH binaries - at separa- 
tions a < ah (equation ^|) - would most likely be found in 
galaxies with low central densities and little gas. Recent 
formation via mergers, and a high ongoing accretion rate 
(assuming that the accreted galaxies also contain BHs), 
would also be propitious. These characteristics constitute 
almost a textbook definition of a cD galaxy, partic ularly a 



centei | which grows as the binary shrinks; at some point 



the binary lies entirely within the hole and its decay ceases 
Using Figure |lj, it is easy to show that the critical sep- 
aration is ~ a few times 10 -4 (assuming that the decay 
stalls when the hole grows to a radius of a few times a). 
Scali ng to M32. this separation is a ga 0.02 pc. and to M87. 



multiple-nucleus cD in a rich galaxy cluster (e.g. Schnci ■ 
der fc Gunn (1982| )). Most cD galaxies are too distant 



for single, much less double, BHs to be detected kincmat- 
ically, although a strong case can be made for dual BHs 
in 3C75, the central radio source in A400. This galaxy ex- 
hibits a pair of radio jets that appear to be emitte d from 
point sourc es separated by ~ 7 kpc in projection (Owen 
l |pc. The gravitational radiation time scale at these et al., 1985 ). But a more likely separation for a binary BH 



separations would be 10 13 yr (M32) and 10 11 yr (M87), 
much longer than the age of the universe. In reality, much 
of the ejected mass in the early stages of the decay comes 
from stars with larger pericenters, allowing the binary to 
shrin k more than this, say by a factor of ~ 2. (This is. 
coinci dentally, roughly the final separation reached in our 



would be the much smaller distance at which the decay is 
expected to stall, 0.01 pc <; a <, 1 pc. Binaries with these 
separations might barely be detectable in nearby galaxies 
by using VLBI techniques to resolve the co mpact radio 
sources associated with the individual BHs (Slee et al. 



1994) 
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Th e existence of uncoalesced BH binaries in galactic nu- 



clei h is a number of interesting consequences. BH ejec- 



similar effect i s seen in the merger simulations of Merritt 
& Cruz (2001). Thus it seems possible that dark matter 



tions would result whenever a third supermassive BH, or 
a second BH binary, is broug ht into the nucleus follow- 
ing a merger ( Valtoncn, 1996| ). In a massive galaxy like 
a cD, a quasi-steady state might be set up in which the 
ejecti on of BHs from the nucleus is matched by the infall 



cusps could be destroyed even in galaxies which manage 
to retain steep cusps in the stellar distribution. 

Could the mechanisms discussed in this paper be rele- 
vant to the low apparent density of dark matter a t the cen- 
ters of dwarf and low-surfacc-brightncss galaxies ( Flores & 



of ne w BHs from the "multiple nuclei" and from previous Primack, 1994|; |dc Blok fc McGaugh, 1997| ; [McGaugh fc de 



ejections. If BH ejections are common, most supermassive 
BHs might be located far from the centers of galaxies, and 



only as p DM 



1HZ 



Blok, 1998t|de Blok et al., 2001| )? The density increases 



the nean mass density of BHs in the universe could be Blok et al., 2001), much flatter than predicted by CDM 



much greater than t he mass density inferred from nuclear 
kincmatical studies (Merritt & Ferrarese, 2001bD . 



at the centers of these galaxies 



6.6. Centers of Dark-Matter Halos 

Cold dark matter (CD M) simulations of the growth of 
structure in the uni verse ( Moore et al., 199§| ; ling, 2000: 
Bullock et al., 2001) predict dark- matter halos with steep 
central density cusps, p ~ r^ 1 , 1 <; 7 <J 2, similar to the 
cusps in our initial models. The dense (baryonic) regions 
in which supermassive BHs first for med were probably lo 



cated near the centers of these halos ( Haehnelt, Natarajan 
fe Recs, 199S ) . Cosmological simulations currently lack 



the resolution to handle co mpact massive objects like BHs, 
but a number of authors (|lpser fc Sikivie 1 1987[ |G ondolc| 



fc Silt, 1999|; jUllio, Zhao fc Kamionkowski, 2001|) have 
investigated the response of pre-existing dark matter ha- 
los to the growth of supermassive BHs. An initially steep 
dark-matter cusp becomes even steeper within the radius 
of influence r ar of the BH, pdm 



2.25 < A < 2.5 



dGondolo fc Silk, 199S| ). This result has been claimed to be 
inconsistent with experimental upper bo unds on annihila- 



tion r adiation from the Galactic center ( Gondolo fc Silk., 
' 1999[ [Bcrtonc, Silk fc Sigl, 2000t |Gondolo, 200q ), implying 



cither that dark matter cusps do not exist, or that current 
ideas about the composition of the dark matter are wrong. 

If mergers of dark-matter halos occurred after the BHs 
were in place, however, the effect of the BHs on the dark 
matter density would be roughly the opposite of what 
these studies assume: the BHs would tend to destroy the 
cusps via ejection of dark matter particles. Just the first 
step in the BH merger process - formation of a hard bi- 
nary via dynamical friction - transfers enough energy to 
the background to convert an r~ 2 cusp into a shallower, 
~ cusp within a radius that contains several times the 
BHs' mass (§||; Figure 0). (This conclusion is independent 
of any uncertainties about Brownian motion, loss-cone re- 
filling, etc., which would be negligible anyway in a dark- 
matter-dominated cusp.) Most of the annihilation radia- 
tion from a putative dark-matter cusp would come from 
a region smaller than this; reducing the cusp slope from 
~ — 2 to ~ — 1 lowers the predicted flux by several orders of 
magnitude ( |Gondolo fc Silk (1999j ), Fig. 2). The fact that 
steep cusps persist in the stellar density in many galaxies, 
including the Milky Way ( Alexander, 1999D , suggests that 
dark matter cusps might also sometimes avoid destruction. 
But the mechanisms discussed above for preserving steep 
cusps in the stellar distribution - e.g. star formation from 
infalling gas - are less applicable to dark matter. Further- 
more we expect the coupling between baryons and dark 
matter to be less than perfect, and stellar cusps would 
themselves inject energy into the dark matter as they spi- 
ralled to the center of the merging dark matter halos; a 



models. It is unlikely that BHs alone are responsible for 
this deficit, however, for several reasons. There is currently 
no evidence for supermassive BHs in these galaxies, and 
scaling relations like the M, — a relation would suggest 
small values of M,. The inferred core radii of the dark 
matter halos are very large, of order 10 2 — 10 3 pc, imply- 
ing an ejected mass that is much greater than any likely 
value of M # . Dwarf and low-surface-brightness galaxies 
are also unlikely to have had active merger histories, at 
least since the era of formation of the BHs. The only pos- 
sibility we can see for destruction of dark-matter cusps on 
the observed scales in these galaxies would be the existence 
of a significant population of condensed objects tha t pre- 



date the dark matter halos, such as primordial BHs (Carr 



1985) 



7. CONCLUSIONS 

1. Mergers of equal-mass stellar systems containing su- 
permassive black holes and steep central density cusps pro- 
duce nuclei with shallow cusps, p ~ r^ 1 , inside of a break 
radius r&, similar to the luminosity profiles observed at 
the centers of bright elliptical galaxies. Most of the evo- 
lution in the central density occurs within a short time, 
~ 10 6 — 10 7 yr, after the black holes form a binary; the 
cusp continues to flatten thereafter as the binary ejects 
stars via the gravitational slingshot. The dependence of 
core properties on black hole mass in observed galaxies is 
shown to be consistent with this formation model. 

2. The merger- induced rotation in the nucleus is reduced 
significantly by the binary as it preferentially ejects stars 
whose angular momenta are aligned with its own. The 
stellar velocity dispersion tensor in the nucleus becomes 
mildly tangentially anisotropic as well, although this effect 
is too small to be easily observed in real galaxies. 

3. Hardening of the black-hole binary takes place effi- 
ciently in our simulations due to the large supply of stars 
provided by the dense cusps, and also to the Brownian 
motion of the binary, which allows it to interact with a 
larger number of stars than if it remained fixed. There 
is no significant dependence of the binary hardening rate 
on number of particles N and the binary's loss cone never 
approaches depletion, in spite of the fact that the decay 
is followed over a factor of ~ 20 in semimajor axis after it 
first forms a hard binary, considerably farther than in ear- 
lier simulations. The hardening rate that we measure, if 
it remained constant, would imply gravitational-radiation 
coalescence in a relatively short time, of order 10 8 years, 
following the merger. 

4. However, we argue that the decay of a black-hole 
binary in a real galactic nucleus would sometimes be ex- 
pected to stall at separations of 0.01 — 1 parsec due to de- 
pletion of the stellar loss cone around a nearly-stationary 
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binary. At these separations, gravitational radiation would 
be ineffective at inducing coalescence and the binary would 
persist indefinitely unless some other physical process 
were able to extract its binding energy. We argue that un- 
coalesced black-hole binaries are most likely to be found 
in the nuclei of cD or other giant elliptical galaxies. 

5. If we artificially combine the two black holes just after 
they form a hard binary, the merger remnant preserves 
its steep, p ~ r -2 density cusp. We propose this as a 
model for the retention of steep cusps in the "power-law" 
galaxies, and suggest that these galaxies experienced their 
last major mergers during the quasar epoch. 

6. Our simulations can also be interpreted as describ- 
ing mergers of dark-matter cusps containing supermassive 
black holes. We argue that the steep cusps predicted by 
cold-dark-matter cosmologies would be destroyed by bi- 
nary black holes in galaxies where the stellar cusps are 
also destroyed. 
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APPENDIX A 



THE COULOMB LOGARITHM 

Estimates of the orbital decay rate due to dynamical friction in were dependent on the Coulomb logarithm In A. 
Here we derive estimates for In A in the two cases of interest: a single massive object near the center of a stellar system 
with a steep density profile; and a sphere of finite size representing a merging cusp. In both cases we find In A » 1 . 

The deceleration of a massive test particle due to dynamical friction is often written in the form (e.g. phandrasckhar 
p94^ ); [Spitier (1987Q ) 

{ ^ )= -*E&M££*m (ad 

II yZ 

where M and v are the mass and the velocity of the test particle, p is the density of light field particles, F(v) is the 
fraction of field particles with velocities less than v, and In A is the Coulomb logarithm that arises in the integration over 
a field of constant density. The force exerted on the test particle by an infinite homogeneous field diverges and a cutoff in 
the form of a maximum impact parameter p max is required. The dependence of (Awi) on p max is logarithmic: 



(A«||) ocln + - In A 

V Pmin 



(A2) 



where p. 

field partic les a re assumed to move along rectilinear orbits (e.g 



is a minimum impact parameter cutoff that also nee ds to be specified in context. In treatm ents where the 

von Neumann (1942|)), the integral 



Chandrasekhar 

leading to (A2) diverges at low impact parameters. The divergence vanishes when the proper Keplcrian trajectories are 



used; field stars of a given relative velocity Vo then transmit a net momentum proportional to In y \ + p 2 nax /pQ with 
Po = GM/Vq, and integration over Vo gives a p m in in equation (A2) of order GM/a 2 with a the ID velocity dispersion 
of the field stars. For instance, when the velocity v of the test star is m uch less than a, p m i n ~ GAI/\/2<j 2 , roughly the 
radius of gravitational influence tq of the massive object (Merritt, 2001). 

A variety of prescriptions can be found in the literature regarding the optimal and most accurate choice for p max , and 
thus for In A. Since all gravitationally-bound structures in the universe have finite extent, the Coulomb logarithm is in 
theory a real p hysi cal quantity subject to calculation if one is ready to abandon several simplifying assumptions that 
enter equation ( |Al| ). In practice, the full-fledged pha se-space integration is cumbersome at best and numerical A^-body 
treatments oft en resort to the fitting of equation (Al) to the dynamical drag observed in simulations. 

Maoz (iyyj) derived an expression for dynamical friction in an inhomogeneous isothermal Maxwellian background which 
rebels 



erf(a) - 1 6(|R-r| - d) 



(A3) 



where a = x ■ (R — r)/|R — r| and x = v/v2cr, while R is the position of the test particle and 0(y) = 1 when y > and 
is zero otherwise. The O-function serves to exclude a finite small volume of radius d around the particle from integration, 
necessary since in Maoz's treatment the field-star trajectories are assumed to be straight lines. 

If the stellar system is spherical and centered on the test particle, the radial and the angular integrals in equation ( |A3| ) 
can be separated 



(Avu) = - 



V2G 2 M /2?r 



erf (a) da 



p(r) 



dr 



(A4) 



The first factor in braces can be identified with the velocity factor F(v) appearing in equation (|A1 



ae a eii(a)da = crf(ir) — xcri '(x) = F(v) 



(A5) 



The second factor in braces encapsulates the dependence of dynamical friction on the radial distribution of field particles. 
The formula becomes 

v 2 J d r 

which can be compared with equation ( |Al[ ) to arrive at a definition of the Coulomb logarithm in terms of an arbitrarily 
chosen fiducial density p: 

p lnA= / ?^-dr. (A7) 



Clearly, Maoz's formula reduces to equation ( A2) if the density is constant in an annulus with inner and outer radii d and 
Pmax and vanishes outside. 

Real stellar systems do not ha ve la rge-radius density cutoffs but the density typically decays as a power law p ~ r~ x . 
The spatial integral in equation (A3) converges for any A > 0, and one is free to take the limit p max — > oo. To illustrate 
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this, we calculate the contribution due to the far field (r > d) particles obeying the simplest power-law profile centered 
on the test particle: 



?W=Po(~) ■ (A8) 



Assume for the moment that p(r) = when r < d 7 reflecting a "hole" in the density cusp due to, e.g., loss-cone depletion 
around a compact massive object. In this case the Coulomb logarithm must be chosen as 

, . f co /r\-*dr 1 p(d) 

leading to dynamical friction that is proportional to the density at the hole's inner edge and inversely proportional to the 
logarithmic slope: 

4.G>MF(v)p {d ) _ 

This result helps expose the inadequacy of conclusions drawn in the context of homogeneous backgrounds where greatest 
contribution to the drag force comes from distant encounters (r ^> d). In our case, the fractional contribution from 
near-field particles at distances d < r < 2d amounts to 1 — 2~ A which is more than 50% when A > 1. Also, while in view 
of the traditional choice In A = p max / 'p m in one is prone to expect p m in "C Pmax, we find that Pmax/Pmin ~ e 1 ^ ~ 1 and 
any meaningful choice for the effective p m ax would reflect neither Chandrasekhar's "average distance between stars" nor 
the "size of the system." 



In the absence of a hole of radius d around the test particle, the integral in (A7) may still diverge. This divergence is 
an artefact of an approximation employed by Maoz whereby stars move along straight lines; it vanishes if exact Kepl erian 
trajectories are co mput ed for the field stars, in which case the effective d is roughly tq — GM/a 2 (Spitzer, 1987). In 



applying equation ( [A 10] ) to the orbital decay of a single BH near the center of a stellar system in §|| we take In A = 1/2 
corresponding to 7 = 2 and use for p(d) the mean density within a radius of 0.01, roughly the radius of gravitational 
influence of the BH. 

The second context in which we applied the dynamical friction formula in §0 was the orbital decay of two finite-density 
spheres representing the original cusps of the merging stellar systems. We model such a test object with a spherical mass 
distribution M(r) and assume that it is much denser than the field environment. According to Gauss' theorem, field 
particles coming to within pericenter distance p from the center of the test object do not inte rac t with the entire object 



but only with a portion of mass M(p). This suggests an immediate modification of equation (A3) 



{Av]l) # T m*-A)p{r)« [ e ^ erf(a) - l] 6 (|R_ r| - d) . (All) 

As an application consider estimating the orbital decay rate for a pair of overlapping Jaffe model galaxies in the final 
stages of a merger proceeding along a circular orbit. Each galaxy has density 



M 



r 



^-ssiUJ + (A12) 

and thus M(r) = Mr/r® and a ~ (GM /2ro) 1 / 2 : . Speed of the test body relative to the other body is twice the circular 
velocity v = 2v c = 2(2- 1 ' 2 a) 1 while a can be expressed terms of the angle 9 between r and R, 

xr sin / , •, 

a =-w^\- (A13) 

Inspection reveals that now, with the test body spread out in space, the limit d — > can be taken inside the integral, 
hence the 0-function is identical to u nity. 



With these substitutions equation (All) can be integrated numerically. We emphasize that treating the test galaxy 
as a rigid body is a crude approximation; in reality there will be an outer tidal radius beyond which the galaxies are 
indistinguishable. The integral depends strongly on the choice of tidal radius. A natural choice is the separation between 
the centers R = a, which yields 

(At,,,) « . (A14) 
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APPENDIX B 
COMPUTING DENSITY PROFILES 

Here we present the algorithms which we used to derive smooth estimates, z/(r) and E(i?), of the particle number 
density and surface density profiles from th e TV-body positions. 

The routines in MAPEL (Merritt, 1994) allow one to derive maximally unbiased estimates of v and E using penalty 
functions that embody the approximate power-law nature of these functions. However the MAPEL routines arc relatively 
slow, and this fact presented difficulties when constructing estimates using the N ~ 10 6 particle data sets consisting of 
superposed iV-body output at several time steps. Kernel based algorithms are faster but potentially more biased; however 



we found them to be adequate f or all but the most steeply rising [y 



! ) profiles and so adopted them here. 



Our derivation follows that in [Merritt fc Tremblay (1994 ). In the absence of any symmetries in the particle distribution, 
a valid estimate of the number density v corresponding to a set of particle positions is 



N 



h 



(Bl) 



where h is the window width and if is a normalized kernel, e.g. the Gaussian kernel: 

1 



K{y) 



(2tt) 3 / 2 



-3/72 



(B2) 



Now imagine that each particle is smeared uniformly around the surface of the sphere whose radius is r, ; typically this 
sphere will be centered on the single or binary BH. If the density profile is actually spherically symmetric, this smearing 
will leave the density unchanged; if not, it will produce a spherically symmetric approximation to the true profile. The 
spherically-symmetrized density estimate is 



z>(r) 



N 1 1 

T-— 



d6 sintf K I - 
h 



d 2 = \r-n 



2 i 2 



2rri cos 9 

where 9 is defined (arbitrarily) from the r,;-axis. This may be written in terms of the angle-averaged kernel K: 



(B3a) 

(B3b) 
(B3c) 



N 



v{r) 



K(r, n ,h) = — 

1 

2 



d9 sin (9 K /i" 1 



2rr,- cos ( 



dn K (h 



Substituting for the Gaussian kernel, we find 

K(r, n,h) 



frr. 



(2tt)3/2 V h 2 



e -(rf + r 2 )/2h 2 sinh(rTi/ ^ 



A better form for numerical computation is 

K(r, n,h) 



/rn 



2(2tt) 3 / 2 \h 2 



-(r,-r) 2 /2h 2 



-(ri+r) 2 /2h 2 



(B4a) 
(B4b) 
(B4c) 

(B5) 

(B6) 



We want to vary the window width with position in such a way that the bias-to- variance ratio of the estimate is relatively 
constant. Let hi be the window width associated with the ith particle. The density estimate based on a variable window 
width is 



N 



(B7) 



The optimal way to vary hi is according to Abramson's (1982) rule: 

hi oc v- a (n), a = 1/2. 



(B8) 
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Since we don't know v(rj) a priori, we compute a pilot estimate of v using a fixed kernel and adjust the hi based on this 
estimate (Silverman, 1986). 

The surface density profile could be computed via simple projection of v(r). Instead, we computed directly from 

the coordinates projected along one axis. The two-dimensional kernel estimate of S(R) in the absence of any symmetries 
is 



i=l 

where K 1 is the two-dimensional Gaussian kernel, 



- |R- Ri 



K\y) = ^e-y 2 '\ 

Z7T 



Now smear each particle uniformly in angle 4> at fixed Ri. The density estimate becomes 

" i i 



d 2 = R 2 + R 2 - 2RR, cos ( 



In terms of the angle-averaged kernel K'\ 



N 



i=l 
1 



h 2 ' 



K'(R, Ri, h) = ^-J K'lhT 1 y/R? + R 2 -2RRiCos < 
= l-e-rt+fWloiRRi/h 2 ) 

Z7T 

where the last expression was derived using the Gaussian kernel; Iq is the modified Bessel function. 



(B9) 

(BIO) 

(Blla) 
(Bllb) 

(B12a) 

(B12b) 
(B12c) 



APPENDIX C 
LOSVD EXTRACTION 

We carried out non-param etric recovery of the LOSVDs N(V) on a dense velocity grid, —Vq < V < Vq, by maximizing 
the penalized log-likelihood (Merritt, 1997) 



log C P [N, Vi\ = > log JV(Vi) -aP[N]-n N(V)dV 



V 



(CI) 



where n is the number of particles inside an aperture, Vi is the line-of-sight projection of the ith particle's velocity, and 
P\N] is a natural choice for the penalty functional that is large for noisy N but assigns zero penalty to a Gaussian function 
(Silverman, 1986) 



P[N] 



Vo 



-Vo 



d 

~dV 



logTV(F) 



dV. 



(C2) 



When a is very large, maximization of log Cp yields a pure Gaussian distribution. For the purpose of extracting Gauss- 
Hermite (GH) moments, we chose a smaller than necessary for smooth LOSVDs, thereby ensuring that non-Gaussian 
substructure of the distributions is not compromised by smoothing. 

Once an LOSVD is available, parameters in the GH expansion can be readily calculated. Define the GH moments of 
N as 

pOO 

(C3) 



N(V)g(w)Hi{w)dV 

-oo 



where Hi are the Hermite polynomials as defined by Gerhard (1993) and the weight function 



g(w) 



-w 2 /2 



27T7o 



V-Vo 

CO 



(C4) 



has three free parameters (70, Vq, ctq). Following van der Marel & Franx (1993), we choose these parameters such that 
ho[N] = 1 and hi[N] = fi2[N] = 0. This can be achieved by minimizing the sum (ho — l) 2 + h\ + h\ as a function of 
(70 j Vo, Co). Once these para meters are uniquely determined, the higher-order GH moments, including h% and h±, can be 
evaluated from equation QC3|) by numerical integration. 
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